Issue
The FFT code below did not give the result similar to scipy
library of Python. But I don't know what's wrong in this code.
import numpy as np
import matplotlib.pyplot as plt
#from scipy.fftpack import fft
def omega(p, q):
return np.exp((-2j * np.pi * p) / q)
def fft(x):
N = len(x)
if N <= 1: return x
even = fft(x[0::2])
odd = fft(x[1::2])
combined = [0] * N
for k in range(N//2):
combined[k] = even[k] + omega(k,N) * odd[k]
combined[k + N//2] = even[k] - omega(k,N) * odd[k]
return combined
N = 600
T = 1.0 / 800.0
x = np.linspace(0, N*T, N)
#y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
y = np.sin(50.0 * 2.0*np.pi*x)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
yf = fft(y)
yfa = 2.0/N * np.abs(yf[0:N//2])
plt.plot(xf, yfa)
plt.show()
This gives:
Solution
All the above comments, i.e. roundoff errors and implementation correctness, are true but you missed an important thing... FFT Cooley and Tukey original algorithm is working only if the number of samples N
is a power of 2. You did notice that
np.allclose(yfa,yfa_sp)
>>> False
for your current input N = 600
, the discrepancies are huge between your output and numpy/scipy. But now, let's use the closest power of two, in this case N = 2**9 = 512
, which gives
np.allclose(yfa,yfa_sp)
>>> True
Wonderful! Outputs are now identical this time, and it can be verified for other powers of 2 (Nyquist criterion apart) sizes of input signal y
. For in depth explanations, you may read accepted answer of this question to understand why numpy/scipy fft functions may allow all N
(with most efficiency when N
is a power of two, and least efficiency when N
is prime) instead of just handling this error as you should have, with something like:
if np.log2(N) % 1 > 0:
raise ValueError('size of input y must be a power of 2')
or even, using bitwise and
operator (a truly elegant test imo):
if N & N-1:
raise ValueError('size of input y must be a power of 2')
As suggested in the comments, if size of the signal could't be modified so easily, zero-padding is definitely the way to go for this kind of sampling issue.
Hope this helps.
Answered By - Yacola
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.