1
$\begingroup$

I am trying to resample a signal using the fourier method in python using np.fft.rfft() and 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:

plot produced by ipython notebook

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..

(see 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:

enter image description here

enter image description here

$\endgroup$
  • 1
    $\begingroup$ If you have S all zeros except for say S[9]=1, what do you get? Alternatively, if S = [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0] does it match? $\endgroup$ – geometrikal Sep 6 '15 at 13:57
  • $\begingroup$ @geometrikal Good question! See edit. $\endgroup$ – CliffordVienna Sep 6 '15 at 14:03
  • 1
    $\begingroup$ How many samples in the resampled signal? Try halving the last value of the dft (corresponds to f = 0.5) before resizing. $\endgroup$ – geometrikal Sep 6 '15 at 14:30
  • $\begingroup$ @geometrikal P*N=128 in all cases (N=16 original samples and P=8 oversampling rate in the resampled signal) $\endgroup$ – CliffordVienna Sep 6 '15 at 14:33
  • 1
    $\begingroup$ Yes I shall write one up $\endgroup$ – geometrikal Sep 6 '15 at 14:39
4
$\begingroup$

If you perform the full DFT of a signal with odd number of samples, the frequency values in the spectrum are a DC coefficient and (N-1)/2 conjugate pairs of sinusoid coefficients.

If you perform the full DFT of a signal with even number of samples, the frequency values in the spectrum are a DC coefficient, (N-1)/2 conjugate pairs of sinusoid coefficients, and a single coefficient corresponding to frequency 0.5.

rfft only keeps the positive frequency part of the spectrum, and the inverse, irfft, presumably recreates the negative part of the spectrum before performing an inverse-DFT.

For example, for an odd size signal with rfft of [0 1 2 3], the full spectrum would be [3 2 1 0 1 2 3].

For an even size signal with rfft of [0 1 2 3 4] it presumably has a flag set somewhere in the code that indicates an even sized signal, and so the full spectrum would be [3 2 1 0 1 2 3 4]. Note the 4 only appears once.

When you pad the rfft result though, e.g. [0 1 2 3 4 0 0 0 0], it probably recreates the spectrum as [0 0 0 4 3 2 1 0 1 2 3 4 0 0 0 0]. That is, the 4 is now repeated twice, and therefore that frequency component is twice as strong.

Therefore, before padding you must halve the last rfft coefficient iff the signal has an even number of samples.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.