Empezamos creando la funciĆ³n plot_spectrum
que toma la parte real positiva de la transformada de Fourier, la normaliza y la dibuja.
Se normaliza el valor absoluto de la amplitud de la seƱal. Es decir, la seƱal con mƔs amplitud tendrƔ una altura de 1 y las demƔs serƔn proporcionales a esta.
import matplotlib
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (12, 5)
matplotlib.rcParams['image.cmap'] = "gray"
import numpy as np
from numpy import sin, linspace, pi
from pylab import plot, show, title, xlabel, ylabel, subplot
from scipy import fft, arange
def plotSpectrum(y, tit=None):
"""
Plots a Single-Sided Amplitude Spectrum of y(t) con un sampling rate de Fs
y(t) contiene la senyal
Fs es el sampling rate, muestras por segundo
taken from http://glowingpython.blogspot.com.es/2011/08/how-to-plot-frequency-spectrum-with.html
"""
n = len(y) # length of the signal
k = arange(n) # preparamos un vector de tantas posiciones como la senyal
T = 1.0 #n/Fs # duraciĆ³n total de la senyal
frq = k/T # two sides frequency range
frq = frq[range(n/2)] # one side frequency range
Y = fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]
plot(frq, abs(Y), 'r') # plotting the spectrum
# titulo si existe
if (tit):
title(tit)
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
def plot_spectrum(y, tiempo, tit=None, orig=False):
"""
Si orig es True, dibujamos la senyal respecto al tiempo
"""
if orig:
subplot(2,1,1)
# dibujamos senyal
plot(tiempo, y)
xlabel('Time')
ylabel('Amplitude')
# dibujamos espectro
subplot(2,1,2)
plotSpectrum(y, tit)
else:
plotSpectrum(y, tit)
show()
def vector_tiempo(Fs=150.0):
""" Devuelve un vector de tiempo equiespaciado cada 1/Fs sobre un segundo
Es decir, Fs es la frequencia de muestreo en muestras/segundo
La senyal siempre durarĆ” un segundo
La senyal siempre tendrĆ” Fs muestras
"""
Ts=1.0/Fs; # periodo de cada muestra
t=arange(0,1,Ts) # vector de tiempo que dura un segundo con Fs muestras
return t
def cuadrada(Fs=200, ff=5):
""" crea una onda cuadrada de frecuencia ff.
Fs es la frecuencia en muestras por segundo
La senyal siempre tendrĆ” una duraciĆ³n de 1s """
zero = np.zeros(10)
zeros = np.zeros(Fs/ff/2)
ones = np.ones(Fs/ff/2)
count = 0
y = []
for i in range(Fs):
if i % Fs/ff/2 == 0:
if count % 2 == 0:
y = np.append(y,zeros)
else:
y = np.append(y,ones)
count += 1
return y
def pulso(Fs=200, d=20):
""" crea un pulso de duracion d (en muestras)
Fs es la frecuencia de muestreo en muestras por segundo
"""
nPulse = d
y = np.ones(nPulse)
y = np.append(y, np.zeros(Fs-nPulse))
return y
def ruido_blanco(Fs):
return np.random.randn(Fs)
Dibujamos una seƱal seno de 5Hz y su espectro. Observamos el $\delta$ en 5Hz.
# frecuencia de la senyal, en Hz
ff = 5;
# la duraciĆ³n siempre es de 1 segundo,
# p.e. si tomamos 500 muestras plot_spectrum
# nos representarĆ” de 0 a 250 Hz
# siempre darle un poco mƔs del doble de la mƔxima frecuencia
t=vector_tiempo(50)
# generamos la senyal, pasamos a rad/s
y = sin(2*pi*ff*t)
# dibujamos senyal y su espectro
plot_spectrum(y, t, tit=u"SeƱal seno 1D a 5 Hz", orig=True)
# onda cuadrada
Fs=100
ff=5
t=vector_tiempo(Fs)
y=cuadrada(Fs, ff)
plot_spectrum(y, t, tit=u"SeƱal cuadrada a 5 Hz", orig=True)
Donde se aprecia la frecuencia fundamental $h_{0}$ a $5Hz$ y el resto de harmĆ³nicos a 15, 25, 35, etc.
# pulso
Fs=200
t=vector_tiempo(Fs)
y=pulso(Fs, 20) # el pulso dura 20 de las 200 muestras
plot_spectrum(y, t, tit=u"Pulso de 10 ms", orig=True)
ObservƔndose una preciosa Sync. Probamos ahora con ruido blanco.
Contiene todas las frecuencias.
Fs=2000
t=vector_tiempo(Fs)
y=ruido_blanco(Fs)
plot_spectrum(y, t, tit=u"Ruido blanco", orig=True)
Probamos ahora con una suma de cosenos para ver si la FFT es capaz de discriminar cada una de ellas.
t=vector_tiempo(100)
# generamos la senyal sumando varias senyales desfasadas
# 5, 18 y 47 Hz, desfasadas 0, pi/3 y pi/17
y = np.cos(2*np.pi*5*t+0) + 2*np.cos(2*np.pi*18*t+ np.pi/3) + np.cos(2*pi*47*t+np.pi/17)
# dibujamos senyal y su espectro
plot_spectrum(y, t, tit=u"Suma de seƱales", orig=True)