Nils, I modified a little bit your code. Try this version, maybe you will find it helpful. I wouldn't bother much with the next power of 2 unless you are going to implement something in real time - FFT calculation with non-power of 2 on a PC is just fine. Also, I would consider periodograms only if the signal changes over time, and you want to be able to track these changes over time. Otherwise, taking the whole signal and doing one FFT may be OK. Take a look at the code below (blue) , which is a modified version of your code (and read the comments). Ivo import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft, fftshift def signal_spectrum(x, Fs): n = len(x) f = np.linspace(-Fs/2, Fs/2, n) X=abs(fftshift(fft(x))) # move the upper spectrum to negative freqs return [f, X] pi = np.pi Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = np.arange(0,L)*T x = 0.7*np.sin(2*pi*50*t)+np.sin(2*pi*120*t) y = x + 2*np.random.randn(len(t)) plt.figure() plt.plot(Fs*t[:50],y[:50]) plt.title('Signal corrupted with zero-mean random noise') plt.xlabel('Time (milliseconds)') NFFT = 2**np.ceil(np.log2(L)) # only if you really want this plt.figure() Y = fft(y,NFFT)/L f = Fs/2*np.linspace(0,1,NFFT/2+1) plt.plot(f,2*abs(Y[:NFFT/2+1])) plt.title('Single-sided amplitude spectrum of y(t)') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|') # Another way of calculating and plotting the spectrum f,Yd = signal_spectrum(y, Fs) plt.figure() # you want to add this to be able to see multiple plots plt.plot(f,Yd) plt.title('Spectrum of the the signal') plt.xlabel('Frequency (Hz)') plt.ylabel('|Y(f)|') plt.grid(True) plt.show() On 26 February 2010 14:05, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal.
I found a Matlab template at
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
Matlab has a function
nextpow2(L)
Is there a similar build-in function in numpy/scipy ?
I tried to convert the m-file into a pythonic form. What is needed to obtain a similar figure of the single-sided amplitude spectrum using numpy/scipy/matplotlib ?
from numpy import sin, linspace, pi from numpy.random import randn from pylab import plot, show, title, xlabel, ylabel from scipy.fft import fft
Fs = 1000. # Sampling frequency T = 1./Fs # Sample time L = 1000 # length of signal t = arange(0,L)*T x = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t) y = x + 2*randn(len(t)) plot(Fs*t[:50],y[:50]) title('Signal corrupted with zero-mean random noise') xlabel('Time (milliseconds)') # # # #NFFT = 2^nextpow2(L); # Next power of 2 from length of y
Y = fft(y,NFFT)/L f = Fs/2*linspace(0,1,NFFT/2+1) plot(f,2*abs(Y[:NFFT/2+1])) title('Single-sided amplitude spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') show()
What can be done in case of nonequispaced data ?
http://dx.doi.org/10.1137/0914081
Thanks in advance
Nils _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user