[SciPy-User] frequency components of a signal buried in a noisy time domain signal

[Ivo Maljevic]

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 at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100226/30758e11/attachment.html>