Tutorial: Transformada de Fourier no MATLAB

Meu ultimo trabalho na disciplina de Análise de Sistemas Lineares foi fazer uma aplicação da transformada de fourier. Usei o matlab para dimensionar e simular um filtro digital para filtrar ruidos em um sinal. Segue abaixo o código:


%% INICIALIZAÇÃO DAS VARIAVEIS E DEFINIÇÃO DO SINAL
fs = 500;               % frequência de amostragem
Ts = 1/fs;              % periodo de amostragem
t = 0:Ts:1-Ts;          % intervalo de amostragem
n = length(t);          % tamanho do vetor
f1 = 6;                 % frequência 1
f2 = 20;                 % frequência 2
A1 = 0.7;               % amplitude 1
A2 = 1;                 % amplitude 2


x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t);
y = awgn(x,2,'measured');
% plotando o sinal
figure(1)
subplot(2,1,1)
plot(t,x),title('Onda senoidal 6Hz + 20Hz'),grid
subplot(2,1,2)
plot(t,y,'r'),title('Sinal com Ruido'),grid

xlabel('Tempo (s)')
ylabel('y(t)')

%% TRANSFORMADA DE FOURIER
Y = fft(y);        % módulo da transformada de fourier
% plota a transformada em troncos
figure(2)
plot(abs(Y));              
grid
xlabel('Frequência?')
ylabel('Amplitude?')
title('Primeira iteração da FFT')

%% NORMALIZANDO E CENTRALIZANDO O ESPECTRO DE FREQUÊNCIA
N = length(x);          % numero de pontos da fft
w = -fs/2:fs/2-1;       % intervalo de frequência centralizado

% a amplitude após fft é uma função do número de pontos,
%technical note 1702

Yshift = fftshift(Y);                      
% plota a transformada normalizada em troncos
figure(3)
plot(w,2*abs(Yshift)/N);
grid
xlabel('Frequência (Hz)')
ylabel('Amplitude')
title('FFT centralizada')

%% CRIANDO UM FILTRO DIGITAL
fil1 = [zeros(1,N/2-f1),ones(1,1),zeros(1,N/2-1+f1)];
fil2 = [zeros(1,N/2+f1),ones(1,1),zeros(1,N/2-1-f1)];
fil3 = [zeros(1,N/2-f2),ones(1,1),zeros(1,N/2-1+f2)];
fil4 = [zeros(1,N/2+f2),ones(1,1),zeros(1,N/2-1-f2)];
fil = fil1+fil2+fil3+fil4;

figure(4)
plot(w,2*abs(Yshift)/N,w,2*fil,'-.r')
grid
xlabel('Frequência (Hz)')
ylabel('Amplitude')
title('Faixa de Corte do Filtro')
axis([-80 80 0 1.1])
legend('Espectro de Frequência','Frequência de Corte')

%% APLICANDO O FILTRO DIGITAL E A TRANSFORMADA INVERSA DE FOURIER
Yf = fil.*Yshift;
yf = ifft(fftshift(Yf));
figure(5)
subplot(2,1,1)
plot(t,real(yf)),title('Sinal após aplicação do filtro'),grid
subplot(2,1,2)
plot(t,x),title('Sinal Original'),grid

xlabel('Tempo (s)')
ylabel('y(t)')

O código é bem auto explicativo, mas ainda acho válido falar algumas coisas:
- fft retorna o espectro de frequência negativo, se todas as componentes do sinal forem reais.
- fftshift(y) desloca o espectro de frequência no tempo. Como o gráfico da fft(y) é simétrico após Fs/2 (Frequência de Nyquist), o ideal é que a representação do espectro de frequência seja centralizado no zero.
- O algorítimo da fft no matlab retorna uma magnitude (Amplitude) de A*N/2, ou seja, a amplitude real da componente de frequência é Anorm = 2*A/N.
- Quando aplicar a ifft, não aplicar sobre o módulo da fft, pois o espectro positivo vai gerar um deslocamento de pi/2 no sinal de saída devido a uma componente complexa. Após a aplicação de ifft de um sinal puramente real, uma componente complexa da ordem de 10^-15 ainda pode estar presente, por isso é importante plotar apenas a parte real.

O filtro foi dimensionado em função das frequências que nós escolhemos, mas, em um caso real, o sinal adquirido teria uma frequência desconhecida, assim o filtro teria que "envolver" os picos mais relevantes do espectro de frequência.

Até a próxima.