UNIVERSIDAD NACIONAL MAYOR DE SAN MARCOS

(Universidad del Per, DECANA DE AMERICA)

FACULTAD DE INGENIERIA ELECTRONICA Y ELECTRICA

E.A.P. INGENIERIA DE TELECOMUNICACIONES
LABORATORIO DE PROCESAMIENTO DIGITAL

LABORATORIO No.7: TRANSFORMADA

FOURIER

DISCRETA

DE

7.1 OBJETIVO
En este Laboratorio, se implementar la Transformada Discreta de Fourier (DFT) para
analizar el espectro de una seal de entrada en tiempo real. Despus de calcular la DFT de un
bloque de muestras de entrada, se calcular la magnitud del resultado para posteriormente
visualizarlo en el computador.
7.2 EQUIPOS Y MATERIALES REQUERIDOS
El equipo de experimentacin necesario para realizar el presente laboratorio es:
1. PC Pentium IV o superior.
2. Sistema Operativo Windows 2000 o XP.
7.3 SOFTWARE REQUERIDO
El siguiente software es necesario para realizar este laboratorio:
1. Matlab V6.0 o Scilab.
2. Libreras y archivos de soporte.
7.4 ASPECTOS TEORICOS
Muchos campos del conocimiento, como Medicina, ptica, Fsica e Ingenieria
Electrnica, utilizan La transformada de Fourier (FT) coo una herramienta de anlisis. En la
prctica, por ejemplo los grupos de trabajo en compresin de datos como el Joint Photographic
Experts Group (JPEG) y los de la Motion Picture Experts Group (MPEG) utilizan una forma
modificada de la FT. Esencialmente nos permite examinar informacin en el dominio de la
frecuencia en vez de hacerlo en el tiempo, que la gente encuentra ms natural para algunos
datos. Por ejemplo, muchos sistemas estreo tienen filas de luces que resplandecen de acuerdo a
la intensidad de la banda de frecuencia. La intensidad de los agudos, por ejemplo, la mayor
intensidad de luz es leda en la fila de luces, creando una barra de luces que suben y bajan de
acuerdo a la msica. Este el tipo de informacin producida por la FT.
La Transformada Discreta de Fourier (Discret Fourier Transform - DFT) es la
transformada en la que nos centraremos, puesto que este trabaja sobre datos discretos.
El mtodo directo de clculo de la DFT de una secuencia finita de N muestras, es a
partir de su definicin:
N 1

X k x nWN kn
k 0

(7.1)

donde: k = 0, 1, 2, , N 1 y WN e

2
N

y la DFT Inversa:

x n

1 N 1
X k WN kn

N k 0

(7.2)

donde: n = 0, 1, 2, .., N 1
Como x[n] puede ser compleja, utilizando (7.1) directamente, se requieren N
multiplicaciones complejas y (N - 1) sumas complejas para calcular cada valor de la DFT. Por lo
tanto, para obtener todos los N valores se requiere en total N2 multiplicaciones complejas y N(N1) sumas complejas.
La mayora de los procedimientos para mejorar la eficiencia de los clculos de la DFT
explotan las propiedades de periodicidad y de simetra de WNkn.
7.5 ACTIVIDADES
7.5.1 DFT DE UNA SECUENCIA SENOIDAL
Para el clculo de la DFT en forma directa, escriba el cdigo de la siguiente funcin
dft_01.m mostrada en el siguiente script:
Almacenar el cdigo de la funcin como:
dft_01.m
function[X]=dft_01(x)
%
% Calculo de la DFT de modo directo.
%
Xsize=length(x);
% Clculo la DFT (Por el camino menos eficiente)
for m=0:Xsize-1
sum=0;
for n=0:Xsize-1
sum=sum+x(n+1)*(cos(2*pi*n*m/Xsize)-i*sin(2*pi*n*m/Xsize));
end
X(m+1)=sum;
end

Utilizando esta funcin, vamos a verificar la respuesta en frecuencia de una secuencia

senoidal utilizando la DFT de una secuencia senoidal, teniendo en cuenta lo siguiente:
Frecuencia de muestreo:
Frecuencia analgica de x(t):
Secuencia de entrada:
Nmero de muestras en el tiempo:
Resolucin de escala de frecuencia:

fs = 1000 Hz.
f = 40 Hz.
x [n] = cos(2**f*n)
N = length(x)
fk = fs/(N-1)

Almacenar el cdigo del programa como:

Lab_07_dft_01.m
% CALCULO DE LA DFT por el mtodo directo
% utilizando la funcin dft_01().m
close all;

%
%

clear all;
clc;
Datos iniciales.
f=40;
% frecuencia analgica en Hz.
fs=1000;
% frecuencia de muestreo en Hz.
Ts=1/fs;
% periodo de muestreo
n=0:Ts:0.50;
Secuencia de prueba
x=cos(2*pi*f*n);
N=length(x);
Calculo de la DFT de x[n]
X=dft_01(x);
Xabs=abs(X);
Xmax=max(abs(X));
Ploteo de las grficas resultantes
figure(1);
plot(x, '.-b');
% Grfica de x[n]en tiempo discreto.
xlabel(sprintf('%6.5f Segundos entre muestra y muestra', Ts));
title('Muestras de x(t)');
Grfica en el dominio de la frecuencia.
figure(2);
plot(abs(X),'.','Color',[0.41,0.26,0.10]);
xlabel(sprintf('La resolucin de frecuencia es de %5.2f Hz entre
muestras', fs/(length(X)-1)));
title('Mdulo de la DFT de x')

La Figura 7.1 presenta las muestras de x(t) expresadas como una

secuencia de muestras con un periodo de muestreo de Ts=0.001 segundo.

Figura 7.1 Muestras de x(t)

Figura 7.2 DFT de x[n]

En la Figura 7.2 se presenta las muestras resultado del clculo de la DFT sobre x[n],
mostrando que el factor de escala entre muestra y muestra ha sido calculado de la siguiente
manera:

fk

fs
N 1

(7.3)

Por otro lado, se observa que la DFT es una secuencia en el tiempo discreto de N puntos
y en el dominio de la frecuencia tambin se obtienen N componentes espectrales equiespaciados
las primeras N/2 componentes representan a las frecuencias positivas y las siguientes N/2
frecuencias a las negativas (parte compleja negativa), demostrando de esta manera que el
espectro resultante es simtrico en la posicin k=N/2.
Ejercicio 7.1
Tomando en cuenta el cdigo utilizado, desarrolle un nuevo cdigo para las mismas
condiciones de seal x(t), periodo de muestreo, y nmero de muestras N, donde muestre una
grafica en el dominio de la frecuencia solo los primeros N/2 puntos de frecuencia, y escala de
frecuencia en Hz.
Seguidamente, considerar ahora la siguiente secuencia y sus condiciones:
x[n] = cos(2*f1*n) + 0.5*cos(2*f2*n)
Frecuencia de muestreo:
Frecuencia analgica de x(t):

fs = 200 Hz.
f1 = 40 Hz. y f2 = 140 Hz.

De acuerdo a los resultados obtenidos, analcelos y confirme si los resultados son los esperados
u observa alguna discrepancia. Justifique su respuesta.
7.5.2 DFT DE UNA SECUENCIA DE DATOS
Desde el punto de vista de las aplicaciones de la DFT, ahora presentaremos la obtencin
de la DFT de una secuencia de datos de una seal muestreada en el tiempo.
En este caso vamos a tener los datos almacenados en un archivo x_a1.dat la misma
que tendremos que incorporar el proceso de carga de datos al programa de la DFT del
experimento anterior utilizando el comando load de Matlab.
Los resultados obtenidos, se observan en la Figura 7.3 que nos presenta los datos en el
tiempo con un periodo de muestreo de Ts=125 mseg. Seguidamente en la Figura 7.4 se muestra
el mdulo de la DFT de los datos cargados.

Figura 7.3 Grfica de los datos x_a1.dat

Figura 7.4 Grfica de la DFT de los datos x_a1.dat

Ejercicio 7.2
Desarrolle un programa que lea el archivo x_a1.dat, tal que pueda verificar lo
mostrado en las Figuras 7.3 y 7.4 e identificar alguna diferencia.
Hallar el rea bajo la curva de la respuesta en frecuencia y discutir que es lo que
representa como informacin de anlisis de la seal en prueba.
Observe el valor de frecuencia de cada uno de los valores mximos que presenta la
DFT. Habr alguna relacin entre los valores de frecuencia de cada uno de estos mximos?
Justifique su respuesta.
7.5.3 EFECTO DE APLICAR VENTANAS
El bloque de la seal de entrada almacenado para el clculo de la DFT genera ciertos
efectos sobre el espectro obtenido debido a los cambios abruptos al inicio y final de la trama de
datos ya que generalmente se procesan seales cuyas frecuencias no son mltiplos enteros de la
frecuencia de muestreo. Para evitar y/o disminuir la aparicin de frecuencias espurias se utiliza
la tcnica de las ventanas, de manera similar como en el caso de los filtros, pero en este caso es
aplicado a los datos obtenidos.
En l siguiente script, se genera una DFT de 512 puntos, donde una funcin ventana,
rectangular (rectwin) se multiplica con las muestras de la secuencia de entrada.
% EFECTO DEL USO DE LAS VENTANAS
close all
clear all
N = 512;
% nmero de puntos de data
M = 512;
% nmero de muestras
Fs = 8000;
% frecuencia de muestreo
f = 200;
% frecuencia en Hz de la seal a analizar
t = [0 : M-1]' / Fs;
% genera un vector de tiempo
x = cos(2*pi*f*t); % genera el tono de 1 KHz
x = [x;zeros(N-M,1)]; % rellena de ceros segn N - M
x = x .* rectwin(N);
% multiplica por la funcin ventana rectangular
s = abs(dft_01(x));
% calcula la modulo de la DFT
s = s / max(s);
% normaliza de modo que el mximo sea 1.0
fp= Fs * [0:N-1] /N;
% genera un vector de frecuencia
figure(1);
plot(x,'.-r')
% grafica de la secuencia x[n]
figure(2);
a=plot(fp(1:N/2),s(1:N/2),'*-k'); % Grafica de la DFT
set(a,'MarkerSize',3)
xlabel('Frecuencia en Hz');
title('Magnitud al cuadrado de FFT');

grid on;
figure(3);
a=plot(fp,s,'*-b');
% Grafica de la DFT en detalle
set(a,'MarkerSize',3)
axis([0 1000 0 1])
xlabel('Frecuencia en Hz');
title('Vista Detallada');
grid on;

El resultado del efecto de la aplicacin de la ventana Rectangular se muestra en las

Figuras 7.5, 7.6 y 7.7.

Figura 7.5 Grfica de la seal cosenoidal de 200 Hz.

Figura 7.6 Magnitud de la DFT para la ventana rectangular.

Figura 7.7 Magnitud de la DFT en detalle

En el listado siguiente se muestra las funciones ventana que Matlab soporta:

Comando
bartlett
barthannwin
blackman
blackmanharris
bohmanwin
chebwin
flattopwin
gausswin
hamming
hann
kaiser
rectwin
tukeywin
triang

Descripcion
- Ventana de Bartlett.
- Ventana modificada de Bartlett-Hanning.
- Ventana de Blackman.
- Ventana de los 4 trminos mnimos de Blackman-Harris.
- Ventana de Bohman.
- Ventana de Chebyshev.
- Ventana Flat Top.
- Ventana gausiana.
- Ventana de Hamming.
- Ventana de Hann.
- Ventana de Kaiser.
- Ventana Rectangular.
- Ventana de Tukey.
- Ventana triangular>

Ejercicio 7.4
Utilizando el cdigo de uso del efecto de las ventanas, ejecute para f=200 Hz y f=250
Hz . Luego analice y discuta los resultados obtenidos.
Seguidamente, cambie la ventana rectwin por una ventana hamming , luego por la
ventana de kiser , y realizar las mismas simulaciones del caso anterior. Cul es el efecto
del uso de la ventana?
7.5.4 EFECTO DEL RELLENO CON CEROS (PADDING)
Si se tiene N (=2n) muestras de una seal, el nmero de coeficientes de la DFT que se
obtendrn empleando la funcin dft_01.m dentro del programa, tambin es N en el intervalo
[0,2]. En algunos casos no se cuenta con una cantidad de datos que sea una potencia de 2, por
lo que existen dos alternativas:
- Adicionar o eliminar datos para llegar a la potencia de 2n ms cercana.
- Adicionar ceros a los datos obtenidos hasta alcanzar la potencia de 2n superior ms
cercana.
Eliminar datos tiene la desventaja de disminuir la resolucin de la DFT, debido a que la
cantidad de datos N con la que se va a operar es menor. Mientras que al adicionar ceros a los
datos puede introducir alguna contaminacin de armnicos, pero este problema puede ser
mejorado con la utilizacin de las ventanas. En MatLab, si el vector original (seal muestreada)
es de longitud M<N, la cantidad de ceros con la que debe ser rellenado el vector es de N-M
As si M = N = 400 indica que no se tiene relleno de ceros. Notar que N es el total de puntos y
N-M es la cantidad de ceros con que se rellenar el arreglo, es claro que N debe ser mayor que
M.

Figura 7.8 Grfica de la funcin coseno para 50 Hz. sin relleno de ceros

Figura 7.9 DFT resultado sin adicin de ceros

Figura 7.10 Grfica de la funcin coseno para 50 Hz. con relleno de ceros.

Figura 7.11DFT resultado con relleno de ceros.

Ejercicio 7.5
Desarrolle un nuevo cdigo utilizando el ejecutado en 7.5.3 modificando lo necesario
en el script de modo que para N=512:
(a) Obtener la DFT reduciendo a 200 los puntos utilizados (M=200, N=200), para
f = 200 Hz y f = 250 Hz y ventana rectwin.
(b) Obtener la DFT rellenando con ceros de manera de obtener una DFT de 512 puntos
utilizados para f = 200 Hz y f = 250 Hz, con rectwin. (N = 512, M = 350).
(c) Obtener la DFT rellenando con ceros de manera de obtener una DFT de 400 puntos
utilizados para f = 200 Hz y f = 250 Hz , con ventana hamming. (N = 512, M = 350).

(d). Compare los resultados obtenidos en los ltimos tres puntos anteriores.
(e) Realice pruebas para valores menores de M (3 casos), f = 250 Hz y ventanas
rectwin, kaiser y hamming. Analice y discuta los resultados.
(f)
En el cdigo del presente experimento, para la DFT, realice los cambios
necesarios para medir los tiempos de ejecucin para distintos valores de N (para 1024 y 2048)
y luego ejecute el mismo proceso pero utilizando la funcin de la FFT, comente los resultados.
Tener en cuenta que la cantidad de procesos que se ejecutan en su computador sean los mismos
para cada prueba, ya que esto puede alterar los tiempos.
7.5.5 TRANSFORMADA DISCRETA INVERSA DE FOURIER - IDFT
La IDFT dada por la ecuacin 7.2, puede ser expresado como:

x n

1 N 1

2
2
X k cos
kn jsen
kn

N k 0
N

(7.4)

donde n = 0, 1, 2, .., N-1.

La siguiente funcin idft_01.m en Matlab permite calcular la IDFT, ero no es tan
rpido como la funcin ifft:
Almacenado como:
idft_01.m
function [x] = idft_01(X)
%
% Uso
%
[x] = idft_01(X);
%
% Esta funcion no es muy eficiente.
% pero demuestra como se realiza el calculo de la IDFT
Xsize = length(X);
% Realizar la reconstruccion
for n=0: Xsize-1
for m=0: Xsize-1
arra(n+1, m+1) =
X(m+1)*(cos(2*pi*m*n/Xsize)+i*sin(2*pi*m*n/Xsize));
end
end
for n=0:Xsize-1
summ = 0;
for m=0:Xsize-1
summ = summ + arra(n+1, m+1);
end
% Division de la longitu de la salida
x(n+1) = summ / Xsize
end

El primer par del lazo anidado crea una matriz llamada arra, el cual almacena cada
valor necesario para la reconstruccin. El segundo par del lazo anidado halla la suma para cada
columna, y divide la longitud de los datos.
Existen muchas formas para hacer ms eficiente a esta funcin. Seguidamente se
demostrara las funciones de la DFT e IDFT. Primero, definiremos una senal x con algunos

valores. Luego, llamamos a la funcin dft_01. Luego que los valores de la DFT son
almacenados en la variable X, la IDFT es calculada y almacenada en x_2. Esta seal tendr la
misma informacin que x. La funcin round consigue montar ms all del punto decimal y la
funcin real ignora la parte imaginaria (el cual es cero).
El archivo de prueba es almacenado como:
Lab_07_idft_01.m
% CALCULO DE LA DFT DE x
%
% Datos de entrada x
x=[7 4 3 9 0 1 5 2];
X=dft_01(x);
% Funcin DFT de x
XReal=real(X);
% Valor real de X
XImag=imag(X);
% Valor imaginario de X
% Grafica de la DFT
figure(1);
plot(XReal,'*r');
figure(2);
plot(XImag,'*b');
%CALCULO DE LA IDFT DE X
x_2 = idft_01(X)
% Funcin IDFT de X
x_2Real=real(x_2);
% Valor real de x
x_2Imag=imag(x_2);
% Valor imaginario de x
% Grafica de la IDFT
figure(3);
plot(x_2Real,'.b');
title('Valores Reales de la IDFT de X');
figure(4);
plot(x_2Imag,'.r');
title('Valores Imaginarios de la IDFT de X')
xx=real(round(x_2))

Resultados:
X=
Column 1
31.0000

Column 2
4.1716 - 5.0711i

Column 3
-1.0000 + 6.0000i

Column 4
9.8284 - 9.0711i

Column 5
-1.0000 - 0.0000i

Column 6
9.8284 + 9.0711i

Column 7
-1.0000 - 6.0000i

Column 8
4.1716 + 5.0711i

Columna 1
7.0000 + 0.0000i

Columna 2
4.0000 + 0.0000i

Columna 3
3.0000 - 0.0000i

Columna 4
9.0000 - 0.0000i

Columna 5
0.0000 - 0.0000i

Columna 6
1.0000 + 0.0000i

Columna 7
5.0000 + 0.0000i

Columna 8
2.0000 + 0.0000i

x_2 =

xx =
Columna 1 hasta el 5

Columna 6 hasta el 8

10

Como se ve, el resultado del redondeo de xx es el mismo que la de x, osea igual a los
datos originales, demostrando de esta manera que la transformada directa e inversa trabajan
juntos como era de esperarse.
Ejercicio 7.6
Utilizando Matlab, construya una seal como un arreglo de al menos 10 valores. Luego
hallar la DFT de esta seal, y a continuacin plotear las amplitudes vs el nmero de muestras.
Tambin hallar la IDFT de la seal transformada y calcular la diferencia entre la seal original y
la reconstruida. En cuanto se aproxima la seal reconstruida con la original? Es esto esperado?.

7.6 CUESTIONARIO.
1.- Suponer, que hemos muestreado una senal a fs = 5000 muestras por segundo, y que tambin
hemos tomado 250 muestras. Despus de realizar la DFT, encontramos que los 10 primeros
resultados son como sigue:
X[k] = 10, 0, 0, 2 +j4, 0, 1, 0, 0, 0, 2 j4
Qu nos dice esto acerca de la frecuencias que estn presentes en nuestra seal de entrada?
(Asumir que otros valores para X[k] superiores a k = 125 son cero).

2.- Una secuencia de entrada, x[n], tiene su transformada de Fourier realizada. El resultado es:
X[k] = {3, 2 +j4, 1, 5 j3, 0, 0, 0, 5 + j3, 1, 2 j4}
(a) Hallar las magnitudes y ngulos de fase. Luego plotear los resultados.
(b) Usted deber notar alguna simetra en su respuesta para la primera parte. Qu clase de
simetra espera usted para la magnitud, y por que?

11