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