2

I'm trying to smooth and interpolate some periodic data in python using scipy.fftp. I have managed to take the fft of the data, remove the higher order frequencies above wn (by doing myfft[wn:-wn] = 0) and then reconstruct a "smoothed" version of the data with ifft(myfft). The array created by the ifft has the same number of points as the original data. How can I use that fft to create an array with more points.

x = [i*2*np.pi/360 for i in range(0,360,30)]
data = np.sin(x)
#get fft
myfft = fftp.fft(data)
#kill feqs above wn
myfft[wn:-wn] = 0
#make new series
newdata = fftp.ifft(myfft)

I've also been able to manually recreate the series at the same resolution as demonstrated here Recreating time series data using FFT results without using ifft

but when I tried upping the resolution of the x-values array it didn't give me the right answer either.

Thanks in advance

Niall

Community
  • 1
  • 1
nrob
  • 861
  • 1
  • 8
  • 22

1 Answers1

5

What np.fft.fft returns has the DC component at position 0, followed by all positive frequencies, then the Nyquist frequency (only if the number of elements is even), then the negative frequencies in reverse order. So to add more resolution you could add zeros at both sides of the Nyquist frequency:

import numpy as np
import matplotlib.pyplot as plt

y = np.sin(np.linspace(0, 2*np.pi, 32, endpoint=False))

f = np.fft.fft(y)
n = len(f)
f_ = np.concatenate((f[0:(n+1)//2],
                     np.zeros(n//2),
                     [] if n%2 != 0 else f[(n+1)//2:(n+3)//2],
                     np.zeros(n//2),
                     f[(n+3)//2:]))
y_ = np.fft.ifft(f_)
plt.plot(y, 'ro')
plt.plot(y_, 'bo')
plt.show()

enter image description here

Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Thanks for the response Jaime, The only problem is that changes the amplitude, which I would like to preserve. Is it as simple as, you have doubled the length of the fft so the amplitude is halved? – nrob Jun 10 '13 at 08:24
  • 1
    Hi Jamie. For the record (as you possibly know) to get the amplitude right you just need to sort out the normalisation. The ifft normalises by the length of the fft, which we have made artificially long. So to correct the amplitude you need to multiply by len(f_)/len(f) – nrob Jun 10 '13 at 08:50