Dgfefbfd

Descargar como docx, pdf o txt
Descargar como docx, pdf o txt
Está en la página 1de 35

fft - Fast Fourier transform (Transformada Rpida de Fourier)

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

es la raz N-esima de la unidad.

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)')

Es difcil identificar los componentes de frecuencia al observar la seal original. Al


convertirla al dominio de la frecuencia, la transformada de Fourier discreta de la seal ruidosa,
se encuentra tomando la transformada rpida de Fourier (FFT):
NFFT = 2^nextpow2(L); % Potencia cercana de 2 de la longitud de y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Grafica del espectro e amplitud de un solo lado.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
La principal razn del porque las amplitudes no son exactamente 0,7 y 1 es por el ruido.
Varias ejecuciones de este cdigo (incluyendo reclculo de y) producen diferentes
aproximaciones a 0,7 y 1. La otra razn es que usted tiene una seal de longitud finita. El
aumento de L de 1000 a 10000 en el ejemplo anterior, producen aproximaciones mucho mejor
en promedio

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.

Tipos de datos Soportados


fft soporta entradas de datos de los tipos double y single. Si se llama a la fft con la
sintaxis y = fft(X, ...), la salida y tiene los mismos tipos de datos que la entrada X.

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.

[3]. FFTW (http://www.fftw.org)

[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.

[5]. Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal Processing, Prentice-Hall,


1989, p. 611.
[6]. Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal Processing, Prentice-Hall,
1989, p. 619.

[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 Transformada Rpida de Fourier (FFT)


Introduccin
La FFT en una dimension
La FFT en multiple dimensiones

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.

Cuando se usan algoritmos de FFT, debe de distinguirse entre la longitud de la ventana y la


longitud de la transformada. La longitud de la ventana es la longitud del vector de datos de
entrada. Este es determinado, por ejemplo, por el tamao de un buffer externo. La longitud de
la transformada es la longitud de la salida, el DFT calculado. Un algoritmo de FFT rellena o
trocea la entrada para realizar la deseada longitud de la transformada. La siguiente figura
ilustra las dos longitudes.
El tiempo de ejecucin de un algoritmo FFT depende de la longitud de la transformada. Es
ms rpido cuando la longitud de la transformada es una potencia de dos, y casi tan rpido
cuando la longitud de la transformada tiene nicamente factores primos pequeos. Es
tpicamente ms lenta para longitudes de transformada que son primos o con factores primos
grandes. Estas diferencias de tiempo, sin embargo, son reducidos significativamente por los
modernos algoritmos de FFT tales como los algoritmos usados en MATLAB. Ajustar la longitud
de estas transformadas para una mayor eficiencia es usualmente innecesario en la prctica.

1.1.2 La FFT en una dimensin


Introduccin
Ejemplo: Anlisis espectral bsico
Ejemplo: Anlisis espectral de una llamada de una ballena
Ejemplo: Interpolacin de datos

1.1.2.1 Introduccin

La funcin fft de MATLAB devuelve la DFT y de un vector de entrada x usando un


algoritmo de transformada rpida de Fourier:
y = fft(x);

En este llamado a la fft, la longitud de la ventana m = length(x) y la longitud de la


transformada n = length(y) son los mismos.
La longitud de la transformada es especificado es especificado por un segundo argumento
opcional:
y = fft(x,n);

En esta llamada a la fft, la longitud de la transformada es n. Si la longitud de x es menor


que n, x es rellenado con ceros hacia adelante para incrementar su longitud a n antes de
calcular la DFT. Si la longitud de x es mayor que n, nicamente los primeros n elementos de x
son usados para calcular la transformada.

1.1.2.2 Ejemplo: Anlisis espectral bsico


La FFT les permitir estimar eficientemente las componentes de frecuencia en la data a partir
de un conjunto discreto de valores muestreados a una tasa fija. Las cantidades relevantes en
un anlisis espectral estn listadas en la siguiente tabla. Para datos basados en el espacio,
remplace las referencias al tiempo con referencias al espacio.

Cantidad Descripcin

x Datos muestreados

m = length(x) Longitud de la data (nmero de muestras)

fs Muestras/unidad de tiempo

dt = 1/fs Incremento de tiempo por muestra

t = (0:m-1)/fs Rango del intervalo de tiempo para la data

y = fft(x,n) Transformada de Fourier Discreta (DFT)

abs(y) Amplitud de la DFT

(abs(y).^2)/n Potencia de la DFT


Cantidad Descripcin

fs/n Incremento de frecuencia (resolucin)

f = (0:n-1)*(fs/n) Rango del intervalo de frecuencia

fs/2 Frecuencia de Nyquist

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

nextpow2 halla el exponente de la potencia de 2 ms cercana mayor o igual a la longitud de la


ventana (ceil(log2(m))), y pow2 calcula la potencia. Usando una potencia de dos para la
longitud de la transformada optimiza el algoritmo de la FFT, aunque en la prctica exista
usualmente una pequea diferencia en el tiempo de ejecucin al usar n = m.

c) Visualizar la DFT. Las grficas de abs(y), abs(y).^2, y log(abs(y))son todas comunes.


Una grafica de la potencia versus frecuencia es denominada un periodograma:
plot(f,power)
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf Periodogram}')
La primera mitad del rango e frecuencia (de 0 a la frecuencia de Nyquist fs/2) es suficiente
para identificar las componentes de frecuencia en la data, desde que la segunda mitad es justo
un reflejo de l primera mitad.

d) En muchas aplicaciones es tradicional centrar el periodograma en 0. La funcin fftshift re


arregla la salida de la fft con un desplazamiento circular para producir un periodograma
centrado en cero:
y0 = fftshift(y); % Re arregla los valores de y
f0 = (-n/2:n/2-1)*(fs/n); % Rango de frecuencia centrado en 0
power0 = y0.*conj(y0)/n; % Potencia centrado en 0
plot(f0,power0)
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf 0-Centered Periodogram}')
El re arreglo hace uso de la periodicidad en la definicin de la DFT (Ver Discrete Fourier
Transform (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

Las componentes de frecuencias


son en su mayora ocultas por la
aleatoriedad en la fase en
valores adyacentes. La tendencia
al alza en la grfica se debe a la
funcin unwrap (desenvolver),
que en este caso suma a la fase
2 con ms frecuencia de lo que
resta.

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.

a) El siguiente cdigo lee, grafica, y reproduce la data:


[x,fs] = auread('bluewhale.au');
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')
title('{\bf Blue Whale Call}')
sound(x,fs)
Una A "trill" es seguido por una serie de B "moans". La llamada B es la ms simple y fcil de
analizar.

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:

m = length(bCall); % Window length


n = pow2(nextpow2(m)); % Transform length
y = fft(bCall); % DFT of signal
f = (0:n-1)*(fs/n)/10; % Frequency range
p = y.*conj(y)/n; % Power of the DFT

d) Grafique la primera mitad del periodograma, hasta la frecuencia de Nyquist:

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}')

La llamada B est compuesto de una frecuencia fundamental alrededor de 17 Hz y una


secuencia de armnicos, con el segundo armnico enfatizado.

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.

a) Aqu est la data que apareci en el artculo de Gauss:

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

Gauss deseaba interpolar un polinomio trigonomtrico de la forma:

b) El siguiente cdigo usa la fft para realizar un equivalente clculo de Gauss:


d = fft(dec);
m = length(dec);
M = floor((m+1)/2);
a0 = d(1)/m;
an = 2*real(d(2:M))/m;
a6 = d(M+1)/m;
bn = -2*imag(d(2:M))/m;
c) Graficando lo interpolado con la data:
hold on
x = 0:0.01:360;
n = 1:length(an);
y = a0 + an*cos(2*pi*n'*x/360) ...
+ bn*sin(2*pi*n'*x/360) ...
+ a6*cos(2*pi*6*x/360);
plot(x,y,'Linewidth',2)
legend('Data','DFT Interpolant','Location','NW')

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.

a) Los astrnomos han tabulado este nmero de casi 300 aos.


load sunspot.dat
year=sunspot(:,1);
relNums=sunspot(:,2);
plot(year,relNums)
title('Sunspot Data')

Aqu est una mirada ms de cerca a los primeros 50 aos.

plot(year(1:50),relNums(1:50),'b.-');
La herramienta fundamental de procesamiento de seales es la FFT, o la Transformada Rpida
de Fourier.

b) Para tomar la FFT de la data de las manchas solares escriba lo siguiente.


Y = fft(relNums);
Y(1)=[];

El primer componente de Y, Y(1), es simplemente la suma de los datos, y se puede quitar.

c) Una grfica de la distribucin de los coeficientes de Fourier (dada por Y) en el plano


complejo es bonito, pero difcil de interpretar. Necesitamos una forma ms til de examinar los
datos en Y.
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');
d) La magnitud compleja al cuadrado de Y se llama la potencia, y una grafica de la potencia
frente a la frecuencia es un "periodograma".

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')

e) La escala en ciclos/ao es algo inconveniente. Podemos graficar en aos/ciclo y estimar la


longitud de un ciclo.

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.

Uso de la FFT para obtener simple graficas para anlisis espectral

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).

% Calculando el nmero de puntos nicos


NumUniquePts = ceil((nfft+1)/2);
% Como FFT es simtrico, deseche la segunda mitad
fftx = fftx(1:NumUniquePts);

d) Despus, calculamos la magnitud de la fft:

% Tomando la magnitud de la fft de x


mx = abs(fftx);

e) Considere el hecho que MATLAB no escala la salida de la fft por la longitud de la entrada:

% Escale la fft por lo que no es una funcin de la longitud de x

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

Ahora, creamos un vector de frecuencias:

% Este es un vector de frecuencias espaciados regularmente


% con un numero de NumUniquePts puntos.
f = (0:NumUniquePts-1)*Fs/nfft;

f) Finalmente, generamos la grfica con un ttulo y etiquetas para los ejes.

% Generando la grafica, titulo y etiquetas.


plot(f,mx);
title('Power Spectrum of a 200Hz Sine Wave');
xlabel('Frequency (Hz)');
ylabel('Power');
g) Juntando todo lo anterior, tendremos el siguiente archivo de MATLAB:

% 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');

La grafica resultante se ver similar a lo siguiente:


El 6.2 agrega un nuevo objeto espectro. Objetos espectro contienen
informacin de los parmetros de un mtodo particular de estimacin espectral (por ejemplo,
. ). Este objeto proporciona una mejor forma para ver y manipular los
parmetros de estimacin espectral. Vea la pgina de referencia del espectro (Spectrum
Objects) y las correspondiente pginas de referencia para el mtodo de estimacin asociado
(Spectral Estimation Method) para ms informacin. El objeto tiene mtodos para evaluar la
densidad espectral de potencia para tcnicas paramtricas y no paramtricas (convencional) y
el espectro cuadrtico medio para las tcnicas no paramtricas. Para las tcnicas de estimacin
espectral subespacio (MUSIC y EIG), sin embargo, el objeto tiene mtodos para calcular la
pseudoespectro. Como un ejemplo para utilizar el objeto del espectro, ver la demostracin
siguiente: Measuring the Power of Deterministic Periodic Signals (Medida de Potencia de
Seales peridicas deterministas).
Prctica 3: Transformada de Fourier en tiempo discreto
DTFT >Transformada de Fourier en Tiempo Discreto

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:

a) La secuencia x[n] puede tener un nmero infinito de puntos.

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.

Al muestrear la DTFT de esta manera se obtiene la expresin correspondiente a la trasformada


discreta de Fourier DFT que en MATLAB se implementa mediante el algoritmo conocido como
FFT (Fast Fourier Transform).

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

N=fix(N); %aproxima a entero redondeando al entero inferior


L=length(x);
if(L>N)
error(' DTFT: numero de muestras, L, debe ser inferior al numero de frec a
calcular N')
end

%
% 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));

En la funcin anterior se realiz un desplazamiento (fftshift) en frecuencias con objeto de que


los resultados de w se den en el intervalo [- , ].

Ejercicios del apartado 1:

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.

DFT y FFT> Algoritmo FFT para calcular la Transformada Discreta de


Fourier

El algoritmo FFT es una manera eficiente de calcular la DFT. En MATLAB la funcin es


X = fft(x,N)

Calcula la FFT de N puntos del vector x.


El resultado X es un vector de nmeros complejos ordenados con ndice k=0,1, ...N-1.
Si no se da el segundo parmetro se considera como N la longitud del vector. Para que el
algoritmo sea eficiente N debe ser potencia de 2.
Si la longitud de x es menor que N, el vector se rellena con ceros. Si es mayor el vector es
truncado.
x = ifft(X)
Calcula la transformada de Fourier inversa del vector X. Tambin se puede especificar el
nmero de puntos N con ifft(X,N)
X = fftshift(x)
Reordena el vector X en orden creciente de frecuencias de tal manera que la componente
continua queda centrada.

Ejercicios del apartado 2:

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.2.2 Cmo estn relacionados los valores de L, N y la resolucin en frecuencias?

3.2.3 Compruebe que sucede en el espectro de la seal si la una secuencia de L= 100


muestras de x[n] se rellena con ceros hasta N=128.

3.2.4(*) Calcule la inversa de la funcin () utilizando la funcin () para recuperar la


seal en el dominio de tiempos.

3.2.5 Suponga que se desea estudiar el contenido en frecuencias usando la FFT, de la


siguiente seal.
x(t) = 0.0472 cos(2)t + 1.5077) + 0.1362 cos(2)t + 1.8769) + 0.4884 cos(2)t -
0.1852) + 0.2942 cos(2)t -1.4488) + 0.1223 cos(2)t).
Cul es su frecuencia fundamental? Qu frecuencia de muestreo debe usarse? Estime un
valor adecuado de N para obtener suficiente precisin en frecuencias. Represente |()| y la
fase de () en funcin de .
3.2.6 Calcule la inversa de la funcin () utilizando la funcin () para recuperar la
seal en el dominio de tiempos.

Enventanado, "Leakage" y resolucin espectral.


Enventanado
Sea la secuencia x[n] = sen(2n/5)

Limitar la secuencia de entrada al intervalo 0, L-l es equivalente a multiplicar la seal de entrada


x[n] por una ventana rectangular w(n) de longitud L= 40 donde

w[n] = l para 0 <= n < L-l


w[n] = 0 para el resto

Ejercicios del apartado 3:

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

w(n) = 1/2(1-cos(2n/(L-l)) para 0<= n <L-l


w(n) = 0 en el resto
3.3.10.- Comparar la dfft de una ventana rectangular y de una ventana de Hamming usando el
mismo valor de L=50

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.

Un sistema como el descrito adolece de dos grandes inconvenientes:

Implica que se conocen las frecuencias presentes en la seal.


La respuesta en frecuencia del filtro pasa baja debe ser cero salvo para w=0.

No obstante tiene su utilidad si se pretende conocer si determinadas frecuencias estn presentes en


la seal y solo interesan la contribucion de estas componentes.

Para seales periodicas de periodo N el FPBJ puede sustiturise simplemente por un acumulador que
suma L puntos

En el caso N=L el comportamiento de este sistema es un filtro pasa baja ideal.

el sistema anterior tiene h[n]= [1, 1, 1, ...1] con L puntos distintos de 0

y su H(w) para L=10 es


Ejercicios del apartado 4:

Considere la seal discreta con periodo N=10.

x[n] = 8 + 10 sen (2/10)n

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.

Qu componentes de frecuencia son distintos de cero? Qu valores toman y por qu?

3.4.3.- Compare los resultados con los obtenidos mediante la fft

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.

Prctica 4: -- Respuesta en frecuencia de sistemas LTI - Muestreo y


Reconstruccin - Filtrado
Sistemas digitales (Filtros digitales)
Sea un sistema LTI descrito por una ecuacin en diferencias lineal de coeficientes constantes.
La salida y[n] en funcin de la entrada x[n], viene dada por

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.

Cuando x[n] = n], se tiene la resuesta al impulso que denominamos h(n).

4.1.1.(*)- Respuesta al impulso de un sistema definido por su ecuacin en diferencias

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.

NOTA La respuesta constante ,H0, cuando n->00, se denomina respuesta estacionaria.

Compruebe grficamente usando la funcin diff (vease help) que la derivada de la respuesta a la
funcin escaln es la respuesta al pulso unitario.

4.1.4.- Respuesta en frecuencia.

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

H(w) = As para |w|< wc

H(w) = 0 para el resto de frecuencias

siendo wc la frecuencia de corte .

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.

La funcin h(t) en el domino de tiempos se obtiene calculando la transformada de Fourier inversa


de H(w) y corresponde a la expresin.

h(t) = (wc*As/ sinc (wc*t/

4.2.1.(*)- En este ejercicio se pretende realizar la reconstruccin de la seal y(t) = cos(2**f0*t)


con f0 =100Hz; seal muestreada con una frecuencia de muestreo fs= 250Hz

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.

Represente la respuesta del filtro h(t).

Considere las muestras y[n] que corresponden al muestreo con Ts = 1/fs

A partir de la secuencia y[n] y la funcin h(t) represente el proceso de reconstruccin que se


corresponde con la aplicacin de un filtro pasa baja con una frecuencia de corte fc=120Hz.

Para comprobar el proceso de reconstruccin se deber representar en un grfico simultneamente la


seal "continua", las muestras tomadas con fs= 250 Hz , y la seal reconstruida.

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

H( w) = constante para w < wo (banda de paso ).


H(w) = 0 para w> wo (banda de rechazo ).

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.3.- A partir de los vectores b y a obtenidos, construya la funcin H(w) siguiendo el


procedimiento del ejercicio 4.1.4 y represente H(w) en funcin de w en el intervalo 0, .
Compruebe que se obtiene el mismo resultado que utilizando la funcin freqz..
4.3.4.(*)- Considere que la seal de entrada est formada for la suma de dos seales sinusoidales de
frecuencia 200 y 800 Hz respectivamente. Usando la funcin filter compruebe que a la salida ha
desaparecido la componente de 800 Hz.
4.3.5.-.- Manteniendo la frecuencia de corte en la funcin ellip modificar uno a uno los
parmetros N, Rp y Rs. Use la funcin freqz para representar las caractersticas H(w) del filtro
y estudie su dependencia con los parmetros anteriores.

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.

REPRESENTACIONES DE FOURIER PARA LAS SEALES


Existen cuatro representaciones distintas de Fourier, cada una aplicable a diferentes tipos de
seales. Estas cuatro clases estn definidas por las propiedades de periodicidad de una seal y si el
tiempo es de tipo continuo o discreto. Las seales peridicas tienen representacin en series de
Fourier. La Serie de Fourier (FS) aplica a seales peridicas de tiempo continuo mientras que la
Serie Discreta de Fourier (DTFS) aplica a seales peridicas de tiempo discreto. Las seales no
peridicas tienen representacin en forma de transformada. Si la seal es continua en el tiempo y no
peridica, la representacin es llamada Transformada de Fourier (FT). Si la seal es discreta en el
tiempo y no peridica entonces la representacin usada es la transformada de Fourier en tiempo
discreto (DTFT). La siguiente tabla ilustra la relacin entre las propiedades de tiempo de una seal
y la representacin de Fourier adecuada.

Tiempo Peridicas No peridicas

Transformada de
Series de Fourier
Fourier
Continuas
( FS )
( FT )

Series discretas de Transformada discreta


Fourier de Fourier
Discretas
( DTFS ) ( DTFT)
La siguiente tabla muestra las relaciones matemticas utilizadas para calcular las representaciones
de Fourier.

Tiempo Peridicas No peridicas


Series de Fourier Transformada de Fourier

Continuas

Series discretas de Fourier Transformada discreta de Fourier

Discretas

La Transformada Discreta de Fourier (DTFS)

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.

Similarmente, dados los coeficientes de una DTFS en un vector llamado X el comando:

>>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.

Considere el siguiente ejemplo:

Determinar los coeficientes DTFS para la siguiente seal:

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 = ones(1,24) + sin( (n * pi / 12) + (3 * pi / 8 ) );

>> X = fft(x)/24;

El resultado terico del ejemplo es el siguiente:

El resultado obtenido mediante los comandos presentados anteriormente es:

X=

Columns 1 through 5

1.0000 0.4619 - 0.1913i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i

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

-0.0000 - 0.0000i -0.0000 - 0.0000i 0 -0.0000 + 0.0000i -0.0000 + 0.0000i

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

-0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.4619 + 0.1913i

Como se puede ver, tres componentes tienen valor diferente de cero.

Un uso comn de la transformada de Fourier, es encontrar las componentes frecuenciales de una


seal en el dominio del tiempo que esta contaminada con ruido. Considrese dos seales senoidales
que tienen frecuencias fundamentales de 50Hz y 120Hz, luego considrese estas seales
contaminadas con ruido aleatorio. Los comandos para generar una seal con las especificaciones
anteriormente mostradas son los siguientes:

>> t = 0:0.001:0.6;

>> x = sin ( 2 * pi * 50 * t ) + sin ( 2 * pi * 120 * t );

>> y = x + 2 * randn ( size ( t ) );

>> plot( 1000 * t (1:50), y (1:50) )

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.

>>Pyy = Y .* conj (Y) / 512;

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))

Para ver todas las muestras y entender la caracterstica de


simetra descrita anteriormente se pueden utilizar los siguientes comandos:

>> f = 1000*(0:511)/512;

>> plot(f,Pyy)

Del espectro de potencia se puede visualizar que las


componentes con mayor frecuencia se encuentran a los 50 y
120 Hz respectivamente. Comprobando as que las seales
de las cuales se formo la seal contaminada con ruido tienen
estas frecuencias fundamentales.

% FFT PARA ANLISIS ESPECTRAL


% Este ejemplo muestra el uso de la funcin FFT para el anlisis espectral.
% Un uso comn de las FFT es encontrar los componentes de la frecuencia de
% una seal enterrada en una seal ruidosa con dominio de tiempo.

% Primero vamos a crear algunos datos. Considerar los datos muestreados


% en 1000 Hz

% Empezaremos formando un eje de tiempo para nuestros datos, que funcione


% desde t=0 hasta t=25 con pasos de 1 milisegundo. Luego formamos una
% seal, x, que contenga ondas seno a 50 Hz y 120 Hz

t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);

% Adicionar algn ruido aleatorio con una desviacin estndar de 2 para


% producir una seal de ruido y. Dar una mirada a esta seal de ruido y
% graficndola.

y = x + 2*rand(size(t));
plot(y(1:50))
title ('Tiempo de Dominio de la Seal Ruidosa')

% Claramente, resulta dificil identificar los componentes de la frecuencia


% de la mirada a esta seal; es por eso que el anlisis espectral es tan
% popular.

% Encontrar la transformada discreta de Fourier de la seal ruidosa y es


% fcil; solo usar la transformada rpida de Fourier (FFT)

Y = fft(y,256);

% Computar la densidad espectral de la energa, una medida de la energa a


% varias frecuencias, usando la compleja conjugada (CONJ). Formar un eje de
% frecuencia para los primeros 127 puntos y usarlo para diagramar el
% resultado. (Los restantes 256 puntos son asimtricos.)

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)')

%____________________________??????

También podría gustarte