Dgfefbfd
Dgfefbfd
Dgfefbfd
Sintaxis
Y = fft(x)
Y = fft(X,n)
Y = fft(X,[],dim)
Y = fft(X,n,dim)
Definiciones
Las funciones Y = fft(x) e y = ifft(X) implementan el par de la transformada y la transformada
inversa sobre los vectores de longitud N mediante:
donde
Descripcin
Y = fft(x) devuelve la Transformada de Fourier Discreta (DFT) del vector x, calculado con
el algoritmo de la Transformada Rpida de Fourier (FFT).
Si la entrada X es una matriz, Y = fft(X) devuelve la transformada de Fourier de cada uno
de las columnas de la matriz.
Si la entrada X es un arreglo multidimensional, fft opera en la primera dimensin no
singleton.
Y = fft(X,n) devuelve la DFT de n-puntos. La fft(X) es equivalente a fft(X,n) donde n
es el tamao de X in the first nonsingleton dimension. Si la longitud de X es menor que n, X es
rellenado con ceros hasta la longitud de n. Si la longitud de X es mayor que n, la secuencia X
es truncada. Cuando X es una matriz, la longitud de las columnas es ajustada en la misma
manera.
Y = fft(X,[],dim) e Y = fft(X,n,dim) aplica la operacin de FFT a travs de la
dimensin dim.
Ejemplo1
Un uso comn de la transformada de Fourier es hallar las componentes de frecuencia de una
seal dentro de una seal ruidosa en el dominio del tiempo. Considere la data muestreada a
1000 Hz. Formemos una seal conteniendo una sinusoide de 50 Hz, de amplitud 0.7 y una
sinusoide de 120 Hz de amplitud 1, corrompida con algn ruido aleatorio de media cero:
Fs = 1000; % Frecuencia de muestreo
T = 1/Fs; % Tiempo de muestreos
L = 1000; % Longitud de la seal
t = (0:L-1)*T; % Vector de tiempo
% Suma de una sinusoide de 50 Hz y una sinusoide de 120 Hz
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoide ms ruido
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
Algoritmos
Las funciones de FFT (fft, fft2, fftn, ifft, ifft2, ifftn) estan basados en una libreria
denominada FFTW [3], [4]. Para calcular la DFT de N-puntos cuando N es compuesta (esto es,
cuando N = N1N2), la libreria FFTW descompone el problema usando el algoritmo de Cooley-
Tukey [1], el cual primero calcula N1 transformadas de tamao N2, y luego calcula N2
transformadas de tamao N1. La descomposicin es aplicada recursivamente tanto para las
DFT de N1 y N2-puntos. Las DFT, hasta que el problema pueda resolverse, utiliza una de las
varias mquinas generados de tamao fijo o "codelets". El codelets a su vez utiliza varios
algoritmos en combinacin, incluyendo una variacin de Cooley-Tukey [5], un algoritmo de factor
prima [6], y un algoritmo de radix dividida [2]. La factorizacin particular de N es elegida
heursticamente.
Cuando N es un nmero primo, la librera FFTW primero descompone un problema de N-
puntos en tres problemas de (N-1)-puntos usando el algoritmo de Rader [7]. Este luego usa la
descomposicin de Cooley-Tukey, descrito anteriormente, para calcular las DFT de (N 1)-
puntos.
Para la mayora de N, las DFT de entrada real requieren aproximadamente la mitad del
tiempo de clculo de la DFT de entrada complejo. Sin embargo, cuando N tiene factores primos
grandes, la diferencia en la velocidad de procesamiento es muy pequea o nada.
El tiempo de ejecucin para la FFT depende de la longitud de la transformacin. Es ms rpido
para las potencias de dos. Es casi tan rpido para las longitudes que tienen nicamente
factores primos pequeas. Es normalmente varias veces ms lento para las longitudes que son
primos o que tienen factores primos grandes.
Note Es posible que pueda aumentarse la velocidad de la fft usando la funcin de utilidad fftw, quien
controla la optimizacin del algoritmo usado para calcular una FFT de una dimensin y tamao
particular.
Referencias
[1]. Cooley, J. W. and J. W. Tukey, "An Algorithm for the Machine Computation of the
Complex Fourier Series,"Mathematics of Computation, Vol. 19, April 1965, pp. 297-301.
[2]. Duhamel, P. and M. Vetterli, "Fast Fourier Transforms: A Tutorial Review and a State of
the Art," Signal Processing, Vol. 19, April 1990, pp. 259-299.
[4]. Frigo, M. and S. G. Johnson, "FFTW: An Adaptive Software Architecture for the
FFT,"Proceedings of the International Conference on Acoustics, Speech, and Signal
Processing, Vol. 3, 1998, pp. 1381-1384.
[7]. Rader, C. M., "Discrete Fourier Transforms when the Number of Data Samples Is Prime,"
Proceedings of the IEEE, Vol. 56, June 1968, pp. 1107-1108.
Ver tambin
dftmtx | fft2 | fftn | fftshift | fftw | filter | freqz | ifft
******************************************************************************
1.1.1 Introduccin
DFT con un milln de puntos es comn en muchas aplicaciones. El procesamiento de
seales e imgenes en aplicaciones modernas serian imposible sin un eficiente mtodo para
calcular la DFT.
La aplicacin directa de la definicin de la DFT (Ver Discrete Fourier Transform (DFT)) ta un
vector de datos de longitud n requiere de n multiplicaciones y n sumas un total de 2n2
operaciones de punto flotante. Esto no incluye la generacin de las potencias de la n-esima raz
compleja de la unidad . Para calcular el DFT de un milln de puntos, un computador capaz de
hacer una multiplicacin y una suma cada microsegundo requerir de un milln de segundos o
alrededor de 11.5 das.
Los algoritmos de Fast Fourier Transform (FFT) tienen una complejidad computacional de
O(nlogn) en lugar de O(n2). Si n es una potencia de 2, una FFT de una dimensin de longitud
n requerir menos que 3nlog2n operaciones de punto-flotante (un nmero de veces de una
constante de proporcionalidad). Para n = 220, es un factor de casi 35,000 veces ms rpido que
2n2.
Las funciones de MATLAB fft, fft2, y fftn (y sus inversas ifft, ifft2, y ifftn, respectivamente) todos
usan algoritmos rpidos de transformada de Fourier para calcular la DFT.
Nota. Los algoritmos de FFT de MATLAB FFT estn basadas en FFTW, "The Fastest Fourier Transform
in the West" (http://www.fftw.org). Ver fft y fftw para mayores detalles.
1.1.2.1 Introduccin
Cantidad Descripcin
x Datos muestreados
fs Muestras/unidad de tiempo
EJEMPLO 2
a) Por ejemplo, genere la siguiente data de x con dos componentes de frecuencia de diferentes
amplitudes y fases contaminados con ruido:
fs = 100; % Frecuencia de muestreo (Hz)
t = 0:1/fs:10-1/fs; % 10 segundos de muestra
x = (1.3)*sin(2*pi*15*t) ... % Componente de 15 Hz
+ (1.7)*sin(2*pi*40*(t-2)) ... % Componente de 40 Hz
+ (2.5)*randn(size(t)); % Ruido Gaussiano;
b) Use la fft para calcular la DFT y la potencia de y:
m = length(x); % Longitud de la ventana
n = pow2(nextpow2(m)); % Longitud de la transformada
y = fft(x,n); % DFT
f = (0:n-1)*(fs/n); % Rango de frecuencia
power = y.*conj(y)/n; % Potencia de la DFT
e) Use las funciones de MATLAB angle y unwrap para crear una grafica de la fase de la DFT:
phase = unwrap(angle(y0));
plot(f0,phase*180/pi)
xlabel('Frequency (Hz)')
ylabel('Phase (Degrees)')
grid on
EJEMPLO 3
1.1.2.3 Ejemplo: Anlisis espectral de la llamada de una ballena
El archive de ejemplo bluewhale.au contiene data de audio de la vocalizacin una ballena
azul del Ocano Pacifico grabada por micrfonos bajo el agua en la costa de California. El
archivo es de la librera de las vocalizaciones de animales mantenido por el Programa de
Investigacin Bioacustica de la Universidad de Cornell.
Nota. La documentacin de los archivos de ejemplo para matemticos de MATLAB estn ubicados en
la subcarpeta \help\techdoc\math\examples de la carpeta raz de MATLAB (matlabroot). Esta
subcarpeta no se encuentra en la ruta de MATLAB en la instalacin. Para usar los archivos de MATLAB
de esta subcarpeta, ya sea que adicione la subcarpeta a la ruta de MATLAB (addpath) o haga la
subcarpeta como la carpeta actual de trabajo (cd).
Debido a que las llamadas de la ballena azul son muy bajas, son raramente audibles a los
humanos. La escala de tiempo en la data esta comprimido por un factor de 10 para elevar el
pitch y hacer la llamada ms claramente audible.
b) Use la grfica anterior para determinar aproximadamente los ndices para el comienzo y el
fin de la primera llamada B. Corrija la base de tiempo por un factor de 10 para acelerar la
data:
bCall = x(2.45e4:3.10e4);
tb = 10*(0:1/fs:(length(bCall)-1)/fs); % Time base
plot(tb,bCall)
xlim([0 tb(end)])
xlabel('Time (seconds)')
ylabel('Amplitude')
title('{\bf Blue Whale B Call}')
c) Use la fft para calcular la DFT de la seal. Corrija el rango de frecuencias por un factor de
10 para acelerar la data:
plot(f(1:floor(n/2)),p(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf Component Frequencies of a Blue Whale B Call}')
EJEMPLO 4
1.1.2.4 Ejemplo: Interpolacin de datos
Este ejemplo demuestra la FFT en un contexto diferente al del anlisis espectral la
estimacin de los coeficientes de un polinomio trigonomtrico que interpola un conjunto de
datos espaciados regularmente. Esta aproximacin de interpolacin de datos es descrito en
[1].
Varias personas descubrieron algoritmos rpidos DFT de forma independiente, y muchas
personas han contribuido a su desarrollo. Un documento de 1965 por John Tukey y John
Cooley [2] es generalmente reconocido como el punto de partida para el uso moderno de la
FFT. Sin embargo, un trabajo de Gauss publicado pstumamente en 1866 [3] (y fechado en
1805) contiene el uso indiscutible de la tcnica de separacin que forma la base de los
modernos algoritmos FFT.
Gauss se interes en el problema de calcular con exactitud las rbitas de los asteroides a
partir de observaciones de sus posiciones. Su artculo contiene 12 puntos de datos sobre la
posicin del asteroide Pallas, a travs del cual desea interpolar un polinomio trigonomtrico
con 12 coeficientes. En lugar de resolver el sistema resultante de 12 por 12 de ecuaciones
lineales a mano, Gauss busc un atajo. Descubri cmo separar las ecuaciones en tres
subproblemas que eran mucho ms fciles de resolver, y luego cmo se recombinan las
soluciones para obtener el resultado deseado. La solucin es equivalente a la estimacin de la
DFT de los datos con un algoritmo FFT.
asc = 0:30:330;
dec = [408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
plot(asc,dec,'ro','Linewidth',2)
xlim([0 360])
xlabel('Ascension (Degrees)')
ylabel('Declination (Minutes)')
title('{\bf Position of the Asteroid Pallas}')
grid on
Referencias
[1]. Briggs, W. and V.E. Henson. The DFT: An Owner's Manual for the Discrete Fourier
Transform. Philadelphia: SIAM, 1995.
[2]. Cooley, J.W. and J.W. Tukey. "An Algorithm for the Machine Calculation of Complex
Fourier Series." Mathematics of Computation. Vol. 19. 1965, pp. 297301.
[3]. Gauss, C. F. "Theoria interpolationis methodo nova tractata." Carl Friedrich Gauss Werke.
Band 3. Gttingen: Kniglichen Gesellschaft der Wissenschaften, 1866.
[4]. Heideman M., D. Johnson, and C. Burrus. "Gauss and the History of the Fast Fourier
Transform." Arch. Hist. Exact Sciences. Vol. 34. 1985, pp. 265277.
[5]. Goldstine, H. H. A History of Numerical Analysis from the 16th through the 19th Century. Berlin:
Springer-Verlag, 1977.
EJEMPLO 5
1.2 Usando la FFT
Esta demostracin usa la funcin FFT para analizar las variaciones en la actividad de las
manchas solares sobre los ltimos 300 aos.
La actividad de las manchas solar es cclica, alcanzando un mximo aproximadamente cada
11 aos. Vamos a confirmarlo. Aqu hay una grafica de una cantidad llamada el nmero
relativa de manchas solares de Zurich, que mide tanto en nmero y tamao de las manchas
solares.
plot(year(1:50),relNums(1:50),'b.-');
La herramienta fundamental de procesamiento de seales es la FFT, o la Transformada Rpida
de Fourier.
n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')
plot(freq(1:40),power(1:40))
xlabel('cycles/year')
f) Ahora graficaremos la potencia frente al periodo por conveniencia (donde periodo=
1/frecuencia). Como esperbamos, existe un muy prominente ciclo con una longitud de
alrededor de 11 aos.
period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');
g) Finalmente, podemos fijar la longitud del ciclo con un poco mas de precisin extrayendo la
frecuencia ms fuerte. Los puntos en rojo ubica este punto.
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;
2 FFT con MATLAB
1) >> X = fft(x)
Hace la FFT del vector x. X es un vector de nmeros complejos ordenados desde k=0...N-1.
Se recomienda que la longitud del vector x sea una potencia de 2. Lo que no se recomienda es
que la longitud de x sea un nmero primo.
2) >> X = fft(x,N)
Especifica el nmero de puntos con el que se quiere hacer la FFT. Si la longitud de x es menor
que N, el vector se rellena con ceros. Si es mayor, el vector es truncado.
3) >> x = ifft(X)
Hace la FFT inversa del vector X.
4) >> x = ifft(X,N)
Especifica el nmero de puntos N con el que quiero hacer la IFFT.
5) >> X = fftshift(X)
Reordena el vector X en orden creciente de frecuencia. Si X es el vector resultante de hacer
una FFT, utilizando esta funcin reordenamos los puntos en funcin de la frecuencia.
2.1 Pregunta
Cmo se puede escalar correctamente la salida de la funcin FFT para obtener una grafica
significativa de la potencia versus frecuencia?
Respuesta
a) Asuma que x es un vector conteniendo su data. Un vector de muestra usado en esta nota
tcnica es una seal sinusoidal de 200 Hz.
% Frecuencia de muestreo
Fs = 1024;
% Vector de teimpo de 1 segundo
t = 0:1/Fs:1;
% Crea una onda senoidal de 200 Hz.
x =sin(2*pi*t*200);
b) Primero, necesita llamar a la funcin FFT. Para una mayor rapidez posible de la s, tendr
que rellenar su datos con ceros suficiente para que su longitud sea una potencia de 2. Al
interior de la funcin FFT lo hace de forma automtica, si se le da un segundo argumento que
especifica la longitud total de la FFT, como se demuestra a continuacin:
% Utilice la prxima potencia de 2 mas alta mayor o
% igual a length(x) para calcular la fft
nfft = 2^(nextpow2(length(x)));
% Tome la fft, rellnelos con ceros hasta
% que length(fftx)sea igual a length(nfft)
fftx = fft(x,nfft);
c) Si nfft es par (que lo ser, si utiliza los dos comandos de arriba), entonces la magnitud de
la FFT ser simtrica, de manera que los primeros (1 + nfft/2) puntos son nicos, y el
resto son simtricamente redundantes. La componente DC de x es fftx(1), y fftx(1 +
nfft/2)es la componente de la frecuencia de Nyquist de x. Si nfft es impar, sin embargo, la
componente de la frecuencia de Nyquist no es evaluado, y el nmero de puntos nicos es
(nfft + 1)/2. Esto puede ser generalizado para ambos casos por ceil((nfft + 1)/2).
e) Considere el hecho que MATLAB no escala la salida de la fft por la longitud de la entrada:
mx = mx/length(x);
% Ahora, tomemos el cuadrado de la magnitud de
% la fft de x el cual ha sido escalado apropiadamente.
mx = mx.^2;
% Desde que dejamos la mitad de la FFT, multipliquemos
& mx por 2 para mantener la misma energa. La componente DC
% y la componente de Nyquist, si este existe, son nicos
% y no debern ser multiplicados por 2.
if rem(nfft, 2) % odd nfft excludes Nyquist point
mx(2:end) = mx(2:end)*2;
else
mx(2:end -1) = mx(2:end -1)*2;
end
% Sampling frequency
Fs = 1024;
% Time vector of 1 second
t = 0:1/Fs:1;
% Create a sine wave of 200 Hz.
x = sin(2*pi*t*200);
% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft= 2^(nextpow2(length(x)));
% Take fft, padding with zeros so that length(fftx) is equal to nfft
fftx = fft(x,nfft);
% Calculate the numberof unique points
NumUniquePts = ceil((nfft+1)/2);
% FFT is symmetric, throw away second half
fftx = fftx(1:NumUniquePts);
% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
mx = abs(fftx)/length(x);
% Take the square of the magnitude of fft of x.
mx = mx.^2;
% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
if rem(nfft, 2) % odd nfft excludes Nyquist point
mx(2:end) = mx(2:end)*2;
else
mx(2:end -1) = mx(2:end -1)*2;
end
% This is an evenly spaced frequency vector with NumUniquePts points.
f = (0:NumUniquePts-1)*Fs/nfft;
% Generamos la grafica, titulo y etiquetas.
plot(f,mx);
title('Power Spectrum of a 200Hz Sine Wave');
xlabel('Frequency (Hz)');
ylabel('Power');
Introduccin.
La transformada de Fourier X(w) de una seal en tiempo discreto x[n] se calcula mediante
la expresin
y su inversa es,
La DTFT X(w) toma valores complejos y es una funcin continua y peridica en w. El periodo
es 2, representndose normalmente en el intervalo [-, ]. Al evaluar numricamente la
DTFT se presentan dos problemas:
b) X(w) es una funcin continua de la frecuencia w y debe ser discretizada para trabajar en
un procesador digital.
Para resolver el primer problema consideraremos que la secuencia de entrada est formada
por un vector de L puntos siendo 0 para los valores comprendidos entre L+ 1 e infinito.
Para el segundo, consideraremos que X(w) se evala en un numero N finito de frecuencias
equidistantes en el intervalo [- , ] con incrementos de 2 /N, es decir se consideran el
conjunto discreto de frecuencias wk = 2k/N con k=0,1,...N-1. Si se elige N lo suficientemente
grande los valores X[2k/N] se aproximan a la funcin X(w) continua origen del muestreo.
Para evitar problemas de muestreo insuficiente se debe elegir N tal que N>L.
Para implementar la dtft usaremos la funcin del el archivo dtft.m que se lista (si es
necesario use help para averiguar cmo funcionan las siguientes lneas)
function [H,W]=dtft(x,N)
% uso: [H W]=dtft(x,N)
% x: muestra de longitud L, se supone que de L+1
% a infinito la muestra toma valor 0.
% N: nmero de frecuencias a evaluar. N debe ser mayor que L.
% H: valores complejos de la DTFT
% W: vector de frecuencias correspondiente a la los valores H calculados
%
% wk=2*pi*k/N con k=0,1,2, ... ,N-1
W=2*pi/N*(0:N-1);
%
medio=ceil(N/2)+1 %aproxima a entero redondeando al entero inferior
%
% evaluamos la DTFT de -pi a pi
%
W(medio:N)=W(medio:N)-2*pi;
W=fftshift(W);
H=fftshift(fft(x,N));
3.1.1 (*) Represente la dtft en mdulo y fase de la seal x[n]= 0.88n*exp ( j(2/5)n), con
L= 40 y N=128.
3.1.2 Compare los resultados y explique qu sucede si se toman valores de N=40, N=64 y
N=1024.
3.1.3 (*) Con N = 128 cambie el valor de L, por ejemplo L=15 y L=128 y comente los
resultados.
3.1.4 Repita los apartados anteriores para la seal x1[n]= exp (j(2 /5)n) y x2[n]=cos((2
/5)n). Explique las diferencias con los apartados anteriores.
3.2.1(*) Sea la secuencia x[n] = cos(0.25n) + cos(0.5n)+ cos(0.52n). Se pide Calcular la DFT
utilizando la funcin matlab (, ) con N=L= longitud de las secuencia x[n] y representar su
mdulo para diferentes valores de nmero de muestras L. Pruebe por ejemplo los siguientes
valores N=16, N=32, N=64, N=128. Indique a partir de qu valor de N son distinguibles las
tres frecuencias de la seal.
3.3.1.(*)-Represente el espectro de las seales w[n], x[n] y del producto y[n]=w[n]*x[n]. Utilice la
funcin fft calculando un numero suficiente de valores ( N=128) para explicar los resultados
anteriores.
3.3.2. Explique, a partir de los espectros anteriores, la relacin del valor mximo obtenido en el eje
de ordenadas al representar |Y(w)| con los parmetros L y N.
Leakage
Una consecuencia del enventanado es que el espectro de la seal no se localiza en una nica
frecuencia. Es decir si tenemos una seal como x(t) = sen wot que solo debera tener una frecuencia
fundamental w=wo, al calcular su FFT tomando una ventana cuadrada, apareceran componentes
adicionales la frecuencia w y su espectro se extiende por todo el intervalo de frecuencias. Este
efecto se conoce como derrame, o bien con el termino ingles "Leakage"
Para comprender este efecto y estimar un valor adecuado para el tamao de la ventana se pide
realizar las siguientes representaciones y estudiar los comportamientos que se presentas en los
siguientes casos.
Ejercicios :
Para todo el ejercicio se considera una seal continua infinita dada por x(t) = sen (2ft) con f =
1KHz. El efecto de aplicar una ventana cuadrada es equivalente a reducir el intervalo de muestreo
en 0 < t < tamao de ventana. Para todos los casos se pide calular x[n], X(w)=DTFT(x[n]), y
X[k]=DFT(x[k]), Representar |X(w)| y |X[k]|.
3.3.3.(*)- Suponga que toma N= 8 muestras considerando el intervalo 0 < t < 1 ms. Cual es la
frecuencia de muestreo Fs?
3.3.4.(*).- Suponga que toma N= 8 muestras considerando el intervalo 0 < t < 0.5 ms. Cual es la
frecuencia de muestreo Fs?
3.3.5.- Suponga que toma N= 24 muestras considerando el intervalo 0 < t < 1.5 ms.
3.3.6.- Suponga que toma N= 64 muestras considerando el intervalo 0 < t < 4 ms.
3.3.7.(*)- Explique que valores de tamao de la ventana son los adecuados para reducir el efecto de
"Leakage"
Resolucin espectral
El enventanado reduce la resolucin espectral (diferencia entre la frecuencia de dos seales para que
pueden ser distinguidas).
Para ello considerar que la seal de entrada viene dada por
3.3.8.- Representar la dtft de esta seal para N=128 y para L=25, 50 y 100. Que relacin hay entre
L y la resolucin en frecuencia.?
Con el fin de reducir el derrame es posible elegir una ventana w(n) cuya dtft W(w) tenga lobulos
laterales ms pequeos, pero esto provoca un aumento en la anchura del lobulo principal, lo que
provoca una disminucin en la resolucin espectral.
3.3.9.- Comprobar este efecto para la seal x[n] anterior usando una ventana de Hamming definida
por
Analizador de Espectros.
Un analizador de espectro es un sistema que permite obtener las freccuencias que estan presentes en
una seal discreta.
El sistema mas simple que nos permite verificar si una seal tiene una componente con frecuencia
wl sera
donde se multiplica la seal de entrada por e-jw1n para desplazar la componente con frecuencia w=w1
al origen de frecuencias w=0 y al aplicar el filtro pasa bajas se obtendr la contribucin de la
componente de la seal con w= w1.
Repitiendo este proceso para cada frecuencia w= wk, se obtendra el espectro del sistema.
Para seales periodicas de periodo N el FPBJ puede sustiturise simplemente por un acumulador que
suma L puntos
3.4.1.- Represente la seal x[n] en un rango adecuado de valores para verificar que es peridica
3.4.2.- Utilizando el sitema descito obentenga los componentes X[k]. Recuerde que el caso
estudiado en el ejemplo corresponde con N=L=10.
3.4.4.- Conocido el espectro exacto de x[n] y del acumulador empleado, explique razonadamente el
funcionamiento del sistema. Por qu el acumulador se puede utilizar como FPBJ?. Represente en
el dominio de frecuencias los diferentes espectros que resultan despus de aplicar el desplazamiento
en frecuencias y su posterior filtrado.
En MATLAB, estas ecuaciones se representan por dos vectores an y bn, expresndose la ecuacin anterior
como
La funcin FILTER(b,a,x) implementa un filtro digital caracterizado por los coeficientes a y b que
filtrar los datos almacenados en x.
Crear los vectores a y b que contengan los coeficientes x[n] e y[n] de la siguiente ecuacin
Genere un impuso unidad con una longitud de 100 elementos, y calcule mediante el comando filter
la respuesta h(n) des sistema anterior. Represente la funcin h[n].
4.1.2.- Resuelva analticamente la ecuacin anterior para n=1,... 5, cuando x[n] = [n] . Suponga las
condiciones iniciales de reposo y[0] = y[-1] = 0;
4.1.3.- Respuesta estacionaria : Calcule mediante la funcin filter la respuesta al impulso h(n)
correspondiente al sistema descrito por la siguiente ecuacin
en el intervalo -10 < n < 100, y la repuesta a la funcin escaln u(n). Represente grficamente la
salida y deternine el valor constante que alcanza la salida (valor que llamaremos H0) para n-> oo.
Verificar, en su caso, que coincide con el valor determinado analticamente.
Compruebe grficamente usando la funcin diff (vease help) que la derivada de la respuesta a la
funcin escaln es la respuesta al pulso unitario.
Para el sistema descrito en el apartado 4.1.3, genere la respuesta a la seal de entrada (impulso en
frecuencia) x[n] = e jn/3 u[n] para -10 < n < 100. Represente grficamente la parte real e imaginaria,
as como la razn entre ambas. Dividendo la respuesta por x(n) se otiene H(w) para w = /3 ,
determine el valor constante, en modulo y fase, de H para n->oo.El valor de H caracteriza en
rgimen estacionario al sistema para una frecuencia . Indica la ganancia o atenuacion y el cambio de
fase para la componente de w = /3 de la seal de entrada
Repita el proceso anterior para distintos valores de w y obtenga una representacion de la respuesta
en frecuencia del sistema H(w).
Compruebe que los valores coinciden con el calculo de H(w) determinado analticamente.
Compare los resultados con los obtenidos mediante la funcin matlab freqz.
Muestreo y Reconstruccin de seales.
La reconstruccin de una seal muestreada a una frecuencia fs cumpliendo el criterio de Nyquist se
realiza de manera ideal con un filtro pasa baja, cuya funcin de transferencia es
Nota: La amplificacin del filtro As permite recuperar sin atenuacin la seal, debido a que el
espectro de la seal original al ser muestreado con un tren de deltas de periodo Ts sufre una
atenuacin que viene dada por un factor 1/Ts.
Represente la seal continua y(t) en el intervalo [-40 ms,40ms] utilizando una variable
independiente t evaluada cada dt= 1/(10*fs) para que conseguir mediante el comando plot (t,y) una
representacin adecuada.
4.2.2. Repita el ejercicio anterior utilizando la funcin conv de Matlab que permite realizar la
convolucin lineal de dos seales.
Se pretende calcular la convolucin y(t) = x(t)*h(t) que se corresponde con aplicar el filtro pasa baja
en el mundo de frecuencia Y(w)=X(w)H(w).
4.2.3. Vari la frecuencia de corte del filtro fc1= 50 Hz, fc2= 100Hz , fc3= 150Hz, fc4=125Hz y
explique los resultados.
Diseo de FILTROS digitales con MATLAB
Los filtros ideales definidos como
no son fsicamente realizables ya que no son causales, sin embargo es posible disear filtros
causales que pueden aproximarse a los ideales con tanta precisin como sea necesaria . Estos filtros
causales presentaran un rizado tanto en la banda de paso como en la banda de rechazo. Por otra
parte la transicin en frecuencia entre ambas bandas no se realizara de forma abrupta. La figura
muestra las caractersticas indicadas
En matlab existen varias funciones que ayudan al diseo de un filtro indicando sus caractersticas.
Una de ellas es ellip (ver help) que permite disear filtros elpticos. Las salidas b y a de esta
funcin son los coeficientes de la ecuacin de diferencias del sistema, definidos en la primera parte
de la prctica. Por otro lado el orden N del filtro se corresponde con el mayor de los valores Na o
Nb empleados en la ecuacin de diferencias.
Unos valores razonables para los parmetros de rizado que deben introducirse en en la funcin ellip
son Rp=0.5 y Rs=20.
4.3.1(*).- Usando la funcin ellip, disear un filtro pasabaja de orden N=4 con frecuencia de corte de 600 Hz
, suponiendo que la seal de entrada ser muestrada a una frecuencia de 8192 Hz..
4.3.2-(*) Mat1ab incorpora la funcin freqz (ver help) que calcula la funcin de transferencia H(w).
4.3.6.-.- Disee un filtro pasa-alta de orden N=6 con frecuencia de corte de 600 Hz, suponiendo que
la seal de entrada ser muestrada a una frecuencia de 8192 Hz. Represente H(w). Considere que la
seal de entrada est formada por la suma de dos seale sinusoidales de frecuencia 200 y 800 Hz
respectivamente. Usando la funcin filter compruebe que a la salida ha desaparecido la componente
de 200Hz,
4.3.7.- Disee un filtro pasa-banda de orden N=6 con frecuencia de paso de 300 y 500 Hz,
suponiendo que la seal de entrada ser muestrada a una frecuencia de 8192 Hz. Represente H(w).
Considere que la seal de entrada est formada por la suma de tres seales sinusoidales de
frecuencia 100,400 y 600 Hz respectivamente. Usando la funcin filter compruebe que a la salida
ha desaparecido la componente de 100 y de 600 Hz.
4.3.8- Disee un filtro suprime-banda de orden N=6 con frecuencia de paso de 300 y 500 Hz,
suponiendo que la seal de entrada ser muestrada a una frecuencia de 8192 Hz. Represente H(w).
Considere que la seal de entrada est formada por la suma de tres seales sinusoidales de
frecuencia 100,400 y 600 Hz respectivamente. Usando la funcin filter se debe comprobar que a la
salida ha desaparecido la componente de 400 Hz.
Transformada de
Series de Fourier
Fourier
Continuas
( FS )
( FT )
Continuas
Discretas
La DTFS es la nica representacin de Fourier que es de valor discreto tanto en el tiempo como en
la frecuencia y de esta manera implcitamente conveniente para una implementacin computacional
en MATLAB. Las expresiones utilizadas para esta representacin son fcilmente implementables en
MATLAB como archivos. Sin embargo los comandos built-in de MATLAB fft y ifft pueden
tambin ser utilizados para evaluar la DTFS. Dado un vector llamado x de longitud N representando
un periodo de una seal peridica x[n]. el comando:
>> X=fft(x)/N
Produce un vector llamado X de longitud N que contiene los coeficientes de la DTFS. Matlab
asume que el periodo evaluado en la seal es desde 0 hasta N-1, de manera que el primer elemento
de x y X corresponden a x[0] y X[0] respectivamente, mientras que los ltimos elementos
corresponden a x[N-1] y X[N-1]. Ntese que la divisin por N es completamente necesaria, debido a
que el comando fft evala la siguiente expresin sin realizar la divisin por N.
>>x=ifft(X)*N
Produce un vector x que representa un periodo de la seal en el tiempo. Ntese que el comando ifft
debe estar multiplicado por N para evaluar la siguiente ecuacin.
Los comandos fft e ifft son computados usando un algoritmo rpido o numricamente eficiente,
conocido como Fast Fourier Transform.
La seal tiene un periodo de 24, de manera que tan solo se hace necesario definir un periodo y
evaluar sobre este periodo la DTFS. Los comandos usados para realizar dicho clculo son:
>> n = 0:23;
>> X = fft(x)/24;
X=
Columns 1 through 5
Columns 6 through 10
-0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i
Columns 11 through 15
Columns 16 through 20
-0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i
Columns 21 through 24
>> t = 0:0.001:0.6;
Es de gran dificultad identificar las componentes de frecuencia mirando la seal original. Sin
embargo al realizar la conversin de esta seal al dominio de la frecuencia, la identificacin de estas
componentes se hace ms sencilla. La conversin de la seal al dominio de la frecuencia se hace
calculando la Transformada R pida de Fourier, tomando para el clculo los primeros 512 puntos de
la seal. El espectro de potencia es una medida de la potencia a varias frecuencias, y este puede ser
calculado con los siguientes comandos.
Para realizar la grfica se puede tener en cuenta que la informacin que aparece en el arreglo Pyy es
por propiedades de la transformada, simtrica con respecto a la frecuencia media, es decir que si
tenemos 512 puntos de muestra, la seal que esta almacenada en el arreglo es simtrica con respecto
a la muestra 256, por lo tanto dibujar las ultimas 256 muestras del arreglo ser completamente
innecesario. De manera que para visualizar el espectro de potencia los comandos deben ser como se
muestran a continuacin:
>> f = 1000*(0:256)/512;
>> plot(f,Pyy(1:257))
>> f = 1000*(0:511)/512;
>> plot(f,Pyy)
t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*rand(size(t));
plot(y(1:50))
title ('Tiempo de Dominio de la Seal Ruidosa')
Y = fft(y,256);
Pyy = Y.*conj(Y)/256;
f = 1000/256*(0:127);
plot(f,Pyy(1:128))
title('Densidad espectral de la energa')
xlabel('Frecuencia(Hz)')
% Enfocar adentro y trazar solamente hasta los 200 Hz. Notar los picos en
% 50 Hz y 120 Hz. Esas son las frecuencias de la seal original.
plot(f(1:50),Pyy(1:50))
title('Densidad espectral de la energa')
xlabel('Frecuencia (Hz)')
%____________________________??????