Skip to main content ↑ ↓

Notice

Please note that most of the software linked on this forum is likely to be safe to use. If you are unsure, feel free to ask in the relevant topics, or send a private message to an administrator or moderator. To help curb the problems of false positives, or in the event that you do find actual malware, you can contribute through the article linked here.
Topic: Why is perfect resampling impossible? (Read 26677 times) previous topic - next topic
0 Members and 1 Guest are viewing this topic.

Why is perfect resampling impossible?

Reply #25 –
Quote
Don't want to drag this thread off topic but you've brought up a question that i've wondered before.

If what you say is correct, would it make more sense upsampling 44.1khz to 88.2khz, instead of upsampling to 48khz when playing back through cards such as Creative's Audigy 2 series.

Thanks
[a href="index.php?act=findpost&pid=237638"][{POST_SNAPBACK}][/a]


Only if the card could actually play the 88.2kHz signal without resampling again, which it won't

Why is perfect resampling impossible?

Reply #26 –
Quote
Only if the card could actually play the 88.2kHz signal without resampling again, which it won't


I was under the impression that Audigy 2 cards resampled everything from 44.1khz downwards but did not resample anything 48khz up to 96khz. I believe I am correct on this  . So in this case it does make sense to go to 88.2khz instead.

Why is perfect resampling impossible?

Reply #27 –
Thank you all very much for your feedback.

Those of you who blame the filtering are probably correct, I suppose. The sampling theorem also requires an ideal lowpass which we do not have. Without an ideal lowpass you can't perfectly reconstruct the original waveform to begin with, so we're stuck.

If you could use an ideal lowpass, the interpolation type should be irrelevant. If I am not mistaken, then ANY form of interpolation only creates additional frequencies above f/2 which could then be removed by the ideal lowpass.

The multiply-everything-with-sinc idea sounds a lot like convolution to me, which would in effect only be another way of going the Fourier route, but I may be completely wrong here.

The Fourier method seems very interesting. I'll see if I can make it work in Matlab. The discrete Fourier transform doesn't offer an optimal lowpass because of discrete frequencies so maybe there will be problems.

Why is perfect resampling impossible?

Reply #28 –
Well, this seems to work pretty good, judging from what I hear and looking at Cool Edit spectrographs:

Code: [Select]
% MATLAB
source = wavread('udial.wav');
spectrum = fft(source);

% length means number of samples
new_length = length(source)*480/441;
padding = round((new_length - length(spectrum))/2);

% ye mägygk:
new_spectrum = fftshift([ zeros(padding,2); fftshift(spectrum); zeros(padding,2) ]);

% there are some non-zero imaginary components which we have to eliminate
resampled = real(ifft(new_spectrum));

wavplay(source,44100)
wavplay(resampled,48000)

wavwrite(resampled,48000,16,'udial_resampled.wav');


The above only works for stereo samples and only does 44.1kHz to 48kHz but should be easily modifiable. For optimal results the input file should have an even number of samples but you could build something that pads in an extra zero sample at the end.

Is udial.wav supposed to have so much high frequency content? I guess that is what trips up bad resamplers. I can hear parts of it and it is damn annoying. I didn't think I could hear that high. My soundcard is a Terratec Aureon Sky (flashed to Prodigy) which shouldn't resample.

Why is perfect resampling impossible?

Reply #29 –
I've been running with the sinc-impulse idea, and it seems that it would be fairly easy to do downsampling with it.

The problem is that you need to filter out everything above the new, lower nyquist limit to avoid aliasing, which the sinc thing does not do. However, if you make the sinc waves wider, you effectively get a low-pass filter.

For example: if you're downsampling by 2, use the function sin(pi x / 2) / (pi x). It is a low-pass filter, and it seems to work pretty well, but I haven't done extensive testing, and I don't know if this would be the "perfect" downsampler. Any thoughts?
"We demand rigidly defined areas of doubt and uncertainty!" - Vroomfondel, H2G2

Why is perfect resampling impossible?

Reply #30 –
Quote
The sampling theorem also requires an ideal lowpass which we do not have.


We have it : perform a DFT, cut the frequencies that you want to remove, then perform the IDFT (Inverse Discrete Fourier Transform). You can also perform a convolution with the needed sinc function, which is exactly the same thing.

Quote
If you could use an ideal lowpass, the interpolation type should be irrelevant. [...] If I am not mistaken, then ANY form of interpolation only creates additional frequencies above f/2 which could then be removed by the ideal lowpass.


I don't think so, look, let's consider this wave file :



If you resample it to 6/5th its original sample rate, you insert one null sample every 5 original samples. Now if you do it with the simplest possible interpolation, which consists in taking the value of the previous sample, you get



It seems clear to me that with some bits of the original wavs put like this, some distortion is introduced below the Nyquist frequency.
An oversampling that consists in doubling the sample rate only introduces frequencies above f/2. But I don't think that it is the case with resampling at other sample rates.
Example, here is the sonogram of udial.wav from 20 to 22050 Hz :



Now, I resample it to 46000 Hz with the lowest quality, and apply an antialias, with SoundForge :



We can see that the interpolation process has created some artifacts well below 22050 Hz. You can remove what you want above with the lowpass that you want, something will always remain below.

Now, if we want to use Fourier Transform in order to resample, we can forget about DFT. All we can do with it is this :



We can apply the brickwall filter, but when we perform the IDFT, we are stuck with one discrete value per original sample. We have still to resample, and we still don't have any perfect algorithm.
Since we need to compute the values taken by the lowpassed function between the samples, we need a complete description of it, not only a discrete one.

We can achieve this with a convolution with the sinc function, that will lead to an ideal lowpass :



This means to sum up the sinc function multiplied with each sample :



Since we can compute sinc(x)=sin(x)/x for any x, we can get the right value for any position.
The only thing we're left with is computing a sum of N values for each sample. N being the total number of samples of the wave file.
We could begin with computing the value of sinc(x) for any x that can be used, but the distances from one original sample to all the new samples can be different for all original samples. So we must calculate the sinc function NxN times.
For a 7 minutes wave file, we have 37,044,000 sample, thus we need to compute 1,372,257,936,000,000 times the value of a given sinc(x). I let math experts search for a possible way to do this.

Why is perfect resampling impossible?

Reply #31 –
Quote
Well, this seems to work pretty good, judging from what I hear and looking at Cool Edit spectrographs:[a href="index.php?act=findpost&pid=237675"][{POST_SNAPBACK}][/a]


If I understand well, you are moving the frequencies inside a new frequency description of the wave, but it seems to me that all your new frequencies will be wrong. If you resample a 1000 Hz sine this way, your padding changes it into a 1044 Hz sine !

Why is perfect resampling impossible?

Reply #32 –
Quote
We could begin with computing the value of sinc(x) for any x that can be used, but the distances from one original sample to all the new samples can be different for all original samples. So we must calculate the sinc function NxN times.
For a 7 minutes wave file, we have 37,044,000 sample, thus we need to compute 1,372,257,936,000,000 times the value of a given sinc(x). I let math experts search for a possible way to do this.
[a href="index.php?act=findpost&pid=237715"][{POST_SNAPBACK}][/a]


With current hardware I think this would take about a year 

Maybe we can have a look at the symbolic result (using the series decomposition of the N*N sinc(x) calculations) and see if something cancels out in the results ?

Why is perfect resampling impossible?

Reply #33 –
The problem is that each of the N values are something like

Sum (sin (xi)/xi, i, 0, N), with xi being a set of N unpredictable x values. I don't see any way to cancel something.

The only hope is that the result becomes unsignificant once |x| gets big enough, because it is majored by 1/|x|. At least it should be possible to estimate the maximum x above which it is unnessesary to go, given the lenght of the wave file, and the value of the lowpass.

Once sum(1/x, x, a, b) gets significantly below the quantization noise (dithering being taken into account), [a,b] being the interval from the outside of the window already calculated until the end of the wav, we can leave it out.

With some luck, maybe we'll only have to calculate 1000 x N sinc values, which is 37,000 times shorter. According to your estimation, it would take 14 minutes.

Why is perfect resampling impossible?

Reply #34 –
Quote
I've been running with the sinc-impulse idea, and it seems that it would be fairly easy to do downsampling with it.
[a href="index.php?act=findpost&pid=237701"][{POST_SNAPBACK}][/a]

FYI, I'm pretty sure that the official (xiph) oggenc uses sinc resampling internally, if that's at all interesting.  However, oggenc's resampling is famously damaging to audible frequencies!  (I assume/hope that's an implementation error...)

Why is perfect resampling impossible?

Reply #36 –
Quote
Quote
I've been running with the sinc-impulse idea, and it seems that it would be fairly easy to do downsampling with it.
[a href="index.php?act=findpost&pid=237701"][{POST_SNAPBACK}][/a]

FYI, I'm pretty sure that the official (xiph) oggenc uses sinc resampling internally, if that's at all interesting.  However, oggenc's resampling is famously damaging to audible frequencies!  (I assume/hope that's an implementation error...)
[a href="index.php?act=findpost&pid=237864"][{POST_SNAPBACK}][/a]


It is an implementation problem.

SRC is for example based on exactly the same system, as are libresample and PPHS.

Why is perfect resampling impossible?

Reply #37 –
Quote
I was under the impression that Audigy 2 cards resampled everything from 44.1khz downwards but did not resample anything 48khz up to 96khz. I believe I am correct on this   . So in this case it does make sense to go to 88.2khz instead.
[a href="index.php?act=findpost&pid=237642"][{POST_SNAPBACK}][/a]


The PCI Audigy 2 variants can only play back multiples of 48 kHz (48,96, possibly 192) without resampling. This can easily be demonstrated (which I have done).

Why is perfect resampling impossible?

Reply #38 –
Quote
If I understand well, you are moving the frequencies inside a new frequency description of the wave, but it seems to me that all your new frequencies will be wrong. If you resample a 1000 Hz sine this way, your padding changes it into a 1044 Hz sine !
[a href="index.php?act=findpost&pid=237716"][{POST_SNAPBACK}][/a]


Well, I just tested it with a single sine at 10kHz and the result was also a sine of 10kHz. A frequency sweep over the entire range from 0 to 22050 Hz produced some artifacts during most of the resampled file at the nyquist frequency of the original 44.1kHz file (gone exactly in the middle).

Originally I had tried myself at explaining what the MATLAB script does, but I couldn't understand what I was saying myself so I deleted it.

Think time and not number of samples. Just for argument's sake: take a 1 second audio snippet, record that at 10Hz. You get 10 samples. DFT those 10 samples and you get 10 fourier coefficients. Now record the same audio snippet at 12Hz -> 12 samples -> 12 DFT coefficients. The coefficients are the spectrum plus a mirrored version. Fftshift helps inject zeros in the middle.

Coefficients 10Hz: a1 a2 a3 a4 a5 a6 a7 a8 a9 a10
Coefficients 12Hz: b1 b2 b3 b4 b5 b6 b7 b8 b9 b010 b11 b12

a1 & b1 are the offset
a2 & b2 represent a sine wave of 1 period
a3 & b3 a sine of 2 periods
...
a6 & b6 a sine wave of 5 periods
a7 = a5 a sine of 4 periods while b7 represents a sine of 6 periods
b8 = b6 a sine of 5 periods
...
a10 a sine of 1 period

b10 a sine of 3 periods
b11 a sine of 2 periods
b12 a sine of 1 period

So if the input and output are both supposed to be 1 second long, then the frequency of an X period sine will be the same --> a2 and b2 represent the same frequency. When you IDFT you get 10 and 12 samples respectively. If both are to represent 1 second of audio, then the first must be assumed to be sampled at 10Hz while the later must be assumed to be sampled at 12Hz.

Since the input doesn't contain frequencies above 5 periods, in the output, all coefficients representing more than 5 periods must be zeroed. --> In the above example b7 must be zero. All other input coeffs must be copied to the output coeffs representing the same number of periods (or frequency).

The above Matlab code should do just that.
- Take input and determine number of samples
- Based on desired resampling, the output, in order to be of the same length, must have N * output_f/input_f samples (the way I built the above code only works for upsampling I think, but you could just as well throw away some of the input coeffs to downsample)
- Construct the spectrum of the output based on the input spectrum and zero all higher frequency components
- Perform IDFT and claim: output is sampled at output_f

edit: I found a little flaw in the resampling process as posted above: The resulting signal will be quieter than the input. This is because from the IDFTs point of view, the total energy in the resampled spectrum is less, but this should be relatively easy to fix.

Another shortcoming is the huge amount of memory needed. Matlab handles audio data as double-precision floats internally. A 30sec 44kHz 16bit stereo audio file takes up roughly 20MB. The complex spectrum twice that much.

Why is perfect resampling impossible?

Reply #39 –
Quote
Take a look at Julius O. Smith III's Digital Audio Resampling Home Page, where ideal resampling is described. The PDF version of the page might be easier to follow than the cluttered web poage.

Laurent de Soras' The Quest For The Perfect Resampler (PDF) describes a similar method with some optimizations and good illustrations of the process.
[a href="index.php?act=findpost&pid=237888"][{POST_SNAPBACK}][/a]

Thank you for those links. They seem very promising.

Why is perfect resampling impossible?

Reply #40 –
Quote
We have it : perform a DFT, cut the frequencies that you want to remove, then perform the IDFT (Inverse Discrete Fourier Transform). You can also perform a convolution with the needed sinc function, which is exactly the same thing.

I don't think this is entirely true: since the DFT produces a discretized spectrum, what happens between one discrete frequency and another?


Quote
If you resample it to 6/5th its original sample rate, you insert one null sample every 5 original samples. Now if you do it with the simplest possible interpolation, which consists in taking the value of the previous sample, you get

I was thinking of it like this:
First step: find a continuous function, that touches all sampled points. An easy way of this would be using linear interpolation between each sample. When you zoom in real close, you will see some zig-zag lines. Now that you have a continuous representation, you can resample by taking new samples at arbitrary points. Most new samples will be along one of thise zig-zag lines. Some coincide with the old samples.

Here's my logic: since all the action is going on in between the old sample points, the only frequencies you can introduce have to be equal or higher than what could previously happen in between 2 sample points. With linear interpolation, you are in effect inserting many sawtooth like signals with period lengths of 1/sampling_rate and additional harmonics above that. Those additional harmonics is what we don't want and have to eliminate with a lowpass. I can see where you're going with the rest of your post, so I may be totally wrong and sawtooths and other interpolations constructed like that have lower frequency components.

Why is perfect resampling impossible?

Reply #41 –
The resampling papers I linked to are a little heavy on the math, so I made a set of graphs to summarize how ideal resampling works. The in:out sampling frequency ratio in this example is 5:6 (5 input samples for every 6 output samples).



The first graph shows the input signal. The next four graphs show the contribution of the first four non-zero input samples to the continuous waveform. Each is a sinc function (sinx/x) scaled by the input sample, with zero-crossings at all input sample points except the middle one. The next graph shows the sum of these four sinc functions, with the original four sample points to show that they match.

The final graph shows the continuous waveform corresponding to the discrete input signal, which is the sum of sinc functions for every input sample. The dots are sample points taken at the output sampling rate. Note how much complexity is hidden in the seemingly innocent input signal (just to be sure I verified that this continuous signal sampled at the input rate matched the input signal).

(I made the diagrams with a little graphic builder I made to experiment with DSP functions)

Why is perfect resampling impossible?

Reply #42 –
Quote
The resampling papers I linked to are a little heavy on the math, so I made a set of graphs to summarize how ideal resampling works. The in:out sampling frequency ratio in this example is 5:6 (5 input samples for every 6 output samples).
[a href="index.php?act=findpost&pid=237942"][{POST_SNAPBACK}][/a]

Thank you very much. This sounds like what others described above.

Quote
Each is a sinc function (sinx/x) scaled by the input sample, with zero-crossings at all input sample points except the middle one.

Does this mean you can get away with a lot less operations?

But I suppose for every target sample, you have to evaluate N times a sinc function and sum it up, since the target sample points will likely not coincide with the zero-crossings (with N being the number of samples in the input).

Why is perfect resampling impossible?

Reply #43 –
Quote
Each is a sinc function (sinx/x) scaled by the input sample, with zero-crossings at all input sample points except the middle one.


Quote
Does this mean you can get away with a lot less operations?

But I suppose for every target sample, you have to evaluate N times a sinc function and sum it up, since the target sample points will likely not coincide with the zero-crossings (with N being the number of samples in the input).


Unless the output rate is a whole multiple of the input rate, you'll only hit zero-crossings occasionally.

The sinc function can be windowed; the shorter the window, the more gradual the rolloff. If you have a sampled sine at just under Nyquist, you'll need a really wide window to properly reconstruct the sine otherwise it'll have amplitude modulation.

I should have mentioned that what I described only works for increasing the sampling rate (where no filtering is necessary). If decreasing the sampling rate, the sinc function should be generated at a lower frequency, i.e. if halving the sampling rate, use sin(0.5x)/(0.5x):



Note how the sinc function used has zero crossings every other input sample (since the output rate is half the input rate). Arbitrary frequency response can be had by filtering the sinc function used.

The input this time is a sine sweep from the the input Nyquist limit to the input rate / 8. You can see above that the output doesn't start to appear until the input sine wave has four samples per cycle, which is when its frequency goes below the Nyquist limit of the output rate.

In both cases the sinc function can be pre-calculated, windowed (since it has infinite length), and sampled at the output interval at several subsample offsets. This reduces processing for each input sample to determining the relative offset to the first output sample that intersects the windowed impulse (which can be done with just an integer multiply then shift), selecting the appropriate set of samples of the sinc function at that offset, and adding each value sequentially to the output buffer.

Why is perfect resampling impossible?

Reply #44 –
Well, in my opinion, programs such as SSRC slow mode or CEP/Audition with pre/post filtering and high quality mode, are able to perfom "virtually" perfect resampling. At last, the difficulty to achieve perfect resampling is just the same difficulty as perfectly brickwall filtering a signal, and these programs filters are very good. Of course, I'm talking of bandlimited interpolation, or what is the same, sinc interpolation, which is the same as "perfect" (upon implementation) interpolation. AFAIK the programs mentioned use this method.

Said that, the only really challenging part of resampling is the filtering part, the rest is not very difficult. AFAIK (I'm not expert in this but I think I know how it works) the easiest way to perform resampling is, summarizing:  first, oversample the signal, then filter everything over original (in case of upsampling) or final (in case of downsampling) fs/2 in order to remove image reflections due to previous process, and then decimate signal.

Say we want to upsample a 32 KHz signal to 48 KHz. The frequency relationship is 1.5, or what is the same, 3/2. The first step, oversampling, is very easy. Just insert 2 zeroes every sample. Now, whe have the signal oversampled 3 times, at 96 KHz. This process leaves intact baseband signal, however it creates image reflections over the original fs/2 (that is, over 16 KHz). The next step is filtering these images, or what is the same, filter everything from 16 KHz and over, using a steep brick wall, linear phase, FIR filter. As a result of this process, signal original amplitude is reduced by a certain factor, that can be easily compensated. Now, we have a signal perfectly upsampled to 96 KHz. The next step, decimation, is achieved just removing 1 every two samples, so at the end of this process whe have the signal decimated by 2, at 48 KHz.

In case of downsampling from, say 48 KHz, to 32 KHz, first we would oversample x2 to 96 KHz, then filter everything over 16 KHz (final fs/2 this time), and then decimate by 3.

For a more general case, the thing is to upsample to the first common multiple of both frequencies, filter, and then decimate.Of course, the trick is to perform the process quickly, perfoming the whole process over chunks of the signal, same way as filtering of signals is usually performed.

Note that this brickwall filtering will, over certain kind of signals, produce time-domain ringing at the frequency cut-off point. In case of upsampling, there will be little or no added ringing, since the original signal must have been already filtered at this same frequency, in order to meet Nyquist requirements for sampling. When downsampling, there may be ringing if there is significant amplitude, impulsive content at the filter frequency cut-off. This ringing can be traded for a gentler roll-off of the brickwall filter, either attenuating the upper part of the original signal spectrum, or allowing some of the image reflections to pass, or a combination of both.

Edit: typos.

Why is perfect resampling impossible?

Reply #45 –
Quote
Take a look at Julius O. Smith III's Digital Audio Resampling Home Page, where ideal resampling is described. The PDF version of the page might be easier to follow than the cluttered web poage.
[a href="index.php?act=findpost&pid=237888"][{POST_SNAPBACK}][/a]


At last! I was looking for that link.

KikeG's idea is simple in terms of understanding the processing, but the above link explains how it's done in practice. Just remember that a sinc function is a low pass filter.


The reason people say that perfect resampling is impossible is this: if you resample (perfectly) to a higher sample rate, and then resample (perfectly) back to the original sample rate, you will not get back to the original data. There's added ringing at the Nyquist limit due to the infinitely long brick wall filter. (You can calculate an infinitely long filter because the input sample is bounded - assuming the signal is zero beyond these bounds means most of this infinite filter can be ignored).

I've often wondered if it's possible to design the filters to ensure that the extra ringing is cancelled out. You can always make the second low pass filter slightly lower than the first, but that's not perfect.

However, it is possible to design resampling filters for integer relationships that do allow perfect up/down-sampling. You ensure that all the original samples (forming every Nth sample in the upsampled version) remain untouched, and the interpolated samples are formed to give the correct result - it is possible (though whether it counts as perfect or not, I haven't checked). Then, to downsample, you simply dump all the interpolated samples and only keep the original ones. Bingo - perfect resampling!


Finally, when upsampling or downsampling where either the source or destination Nyquist frequency is in the audible range, it's unwise to use a perfect brick wall filter - the ringing is easily audible. Try it with a long FFT in Cool Edit on something like Castanets. It sounds horrible. It's better to move the cut off frequency slightly lower, and to use a gentler filter in order to reduce the time domain ringing.

Cheers,
David.

Why is perfect resampling impossible?

Reply #46 –
Just say, again, that when upsampling a real-world signal (not just an artificial delta pulse), this signal will have been already filtered at the ADC at Nyquist frequency (or before), so additional filtering with a very steep filter at this Nyquist frequency will add very little or no ringing, since there will be little content left at this frequency.

In this same case of upsampling a real-world signal, downsampling it then to original frequency, will result in different sample information, but that information in fact will be the same signal, with a little bit of lineal phase displacement, or what is the same, a small time difference, due to the processing performed. If you could perfectly time-align those signals, and substract them, the difference would be very close to zero, in practice probably below noise floor of original signal, even at Nyquist frequency.

However, in the different case of just downsampling a real-world signal, you are removing information from it, and removing frequencies that were present, so there can be ringing for the kind of signals I explained at my previous post, and of course, upsampling the signal to the original frequency you will have lost these frequencies so there will definitely will be a difference.

Why is perfect resampling impossible?

Reply #47 –
Quote from: KikeG on
In this same case of upsampling a real-world signal, downsampling it then to original frequency, will result in different sample information, but that information in fact will be the same signal, with a little bit of lineal phase displacement, or what is the same, a small time difference, due to the processing performed.

Apologies for replying to such an old post, but I don't think that's right. There's no inherent time difference due to processing—providing that the resampling filter is linear-phase, the signals can match perfectly time-wise. I know that many resamplers do introduce sub-sample time-shifts but this is probably the implementation discarding the wrong samples during downsampling; if you discard the right ones, correct phase is maintained.

A major factor here though is quantisation error: to stand a hope of reconstructing the original signal from the intermediate signal, the latter must have higher bit-depth.
 

Why is perfect resampling impossible?

Reply #48 –
Perfect resampling would require a filter of infinite length. You could never succeed in completing the first sample.

As i'm sure somebody's pointed out that you can be "perfect to quantization noise level or below" which is appropriate.
-----
J. D. (jj) Johnston

Why is perfect resampling impossible?

Reply #49 –
Quote from: Woodinville on
As i'm sure somebody's pointed out that you can be "perfect to quantization noise level or below" which is appropriate.
A colleague mastering engineer recently claimed to have measured his Weiss Saracon SRC to be bit-transparent for a 24/44.1->24/96->24/44.1 roundtrip. I was rather surprised since there's low-pass filtering and probably re-dithering involved. Assuming the source material contained some energy close to the Nyquist frequency, wouldn't bit-transparency be quite unlikely ?