I am trying to resample a signal using the fourier method in python using
np.fft.irfft(). (I'm not using
scipy.signal.resample() because in my final application I need to do more to my signal in the frequency domain than what can be accomplished by the
window argument of
scipy.signal.resample().) I have reduced the problems I see to this simple test case (ipython notebook):
%pylab inline N = 16 P = 8 S = [0.49, 1.68, 0.78, 0.05, 0.21, 0.29, 1.34, 1.17, 16.73, 48.78, 16.90, 0.62, 0.40, 1.60, 0.86, 0.57] assert len(S) == N figure(figsize=(14, 5)) plot(np.arange(N), S, '.') dft = np.fft.rfft(S) plot(np.arange(N), np.fft.irfft(dft), '-') dft *= P dft.resize(P*N/2+1) plot(np.arange(P*N)/P, np.fft.irfft(dft), '-') xlim(6.5, 11.5) ylim(-5, 55);
And here is what I get back:
As you can see, the result looks almost correct, but is nevertheless a bit off. (The resampled red curve does not go through the original sample points, even though there are sample points at the same positions in time as in the original signal.)
The greeen line shows that everything is still ok before extending and scaling the
dft array. So the problem seems to be in the actual resampling:
dft *= P dft.resize(P*N/2+1) plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')
I don't see how this result is mathematically possible because other than scaling of time and amplitude axis, adding zeroes to the high-frequency end of the
dft array should only add zero terms to the fourier series expansion of the signal. So the signal should still pass through the original sample values.
But the small error I see does not seem to be related to scaling: sometimes the resampled waveform lags behind the original (t=8) and sometimes it is faster (t=10), sometimes it is larger in magnitude (t=7,9,11) and sometimes it is smaller (t=8,10).
I'm playing with this for hours now and can't find what is wrong with my code.. Maybe I just don't see it because I'm already staring at it for so long now..
np.fft documentation for the definitions used by the python FFT implementations)
Edit (in response to comment): When I replace S with an impulse at t=9, or a fast impulse train, I get the following: