Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Why do scipy and numpy fft plots look different?

I am currently doing some spectrum analysis for a piece of coursework, though we haven't been explicitly taught Fourier transforms yet. I have been playing around with the various fft algorithms in scipy and numpy on some data that I know what the answer should look like

In this case its an AM signal at 8kHz carrier frequency and 1kHz modulated sine wave on top so should have 3 clear peaks on an fft

When applying scipy.fftpack.rfft and numpy.fft.rfft I get the following plots respectively:

Scipy:

enter image description here

Numpy:

enter image description here

While the shape of the 2 FFTs are roughly the same with the correct ratios between the peaks, the numpy one looks much smoother, whereas the scipy one has slightly smaller max peaks, and has much more noise.

I'm assuming this is largely down to different applications of a Discrete Fourier Transform algorithm, and have seen other articles about how the scipy implementation is faster in run time. But I was wandering what it is that specifically causes the difference, and which one is actually more accurate?

EDIT: Code used to generate plots:

data = pd.read_csv("./Waveforms/AM waveform Sine.csv", sep = ',', dtype = float)

data = data.as_matrix()
time = data[:,0]
voltage = data[:,1]/data[:,1].max() # normalise the values

#scipy plot:
plt.figure()
magnitude =  scipy.fftpack.rfft(voltage)
freq = scipy.fftpack.rfftfreq(len(time),np.diff(time)[0])
plt.figure()
plt.plot(freq, np.absolute(magnitude), lw = 1)
plt.ylim(0,2500)
plt.xlim(0,15)

#numpy plot
magnitude = np.fft.rfft(voltage)
freq = np.fft.rfftfreq(len(time),np.diff(time)[0])

plt.figure()
plt.plot(freq, np.absolute(magnitude), lw = 1)
plt.ylim(0,2500)
plt.xlim(0,15)
like image 473
falcoso Avatar asked Nov 15 '17 12:11

falcoso


People also ask

What does Scipy FFT return?

fftpack. fft. Return discrete Fourier transform of real or complex sequence.

What is FFT in Scipy?

Fourier Transforms ( scipy.fft ) Fast Fourier transforms. 1-D discrete Fourier transforms. 2- and N-D discrete Fourier transforms.

What is the difference between FFT and RFFT?

fft returns a 2 dimensional array of shape (number_of_frames, fft_length) containing complex numbers. For np. fft. rfft returns a 2 dimensional array of shape (number_of_frames, ((fft_length/2) + 1)) containing complex numbers.


1 Answers

From NumPy's doc for rfft:

Returns:

out : complex ndarray

The truncated or zero-padded input, transformed along the axis indicated by axis, or the last one if axis is not specified. If n is even, the length of the transformed axis is (n/2)+1. If n is odd, the length is (n+1)/2.

It is not written explicitly but the "transformed data" is here complex.

From SciPy's doc for rfft

z : real ndarray

The returned real array contains:

[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))]              if n is even
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))]   if n is odd

Conclusion: the storage is different.

For a starter, look at the length of magnitude, it will be different in both cases. I give an example below for clarity:

In [33]: data = np.random.random(size=8)

In [34]: np.fft.rfft(data)
Out[34]: 
array([ 3.33822983+0.j        ,  0.15879369+0.48542266j,
        0.00614876+0.03590621j, -0.67376592-0.69793372j,  1.51730861+0.j        ])

In [35]: scipy.fftpack.rfft(data)
Out[35]: 
array([ 3.33822983,  0.15879369,  0.48542266,  0.00614876,  0.03590621,
       -0.67376592, -0.69793372,  1.51730861])

The first element in both cases is the so-called "DC component" (the mean of the signal).

Then, you can recognize in the SciPy version the succession of real and imaginary parts of the NumPy version.

like image 154
Pierre de Buyl Avatar answered Sep 16 '22 18:09

Pierre de Buyl