UNIVERSIDADE FEDERAL DE OURO PRETO
INSTITUTO DE CIÊNCIAS EXATAS E APLICADAS – ICEA
DEPARTAMENTO DE ENGENHARIA ELÉTRICA
Transformada Rápida de Fourier
Disciplina: Processamento digital de sinais
Professor: Fabricio Javier Erazo Costa
Alunos: Elissane Clementina Pinto – 16.1.8093
Laio da Silva Pires – 15.1.5815
Ramon dos Reis Ferreira – 15.1.8393
Introdução
A transformada de Fourier discreta (TFD) desempenha um papel importante na
análise, projeto e implementação de algoritmos e sistemas de processamento em
tempo discreto de sinais, pois as propriedades básicas da transformada de Fourier de
tempo discreto e da transformada de Fourier discreta tornam particularmente
convenientes a análise e o projeto de sistemas no domínio de Fourier.
O foco principal do capítulo é uma classe de algoritmos particularmente eficiente para
a implementação computacional da TFD de N pontos. Coletivamente, esses algoritmos
eficientes são chamados de algoritmos de FFT (do inglês Fast Fourier Transform).
Para alcançar a máxima eficiência, os algoritmos de FFT precisam calcular todos os N
valores da TFD. Quando queremos valores da TFD em apenas algumas frequências
no intervalo , outros algoritmos podem ser mais eficientes e flexíveis,
embora sejam menos eficientes do que os algoritmos de FFT para o cálculo de todos
os valores da TFD.
Existem muitas maneiras de medir o custo e a eficiência de uma implementação ou de
um algoritmo, e uma avaliação final depende tanto da tecnologia de implementação
disponível quanto da aplicação pretendida. Aqui será utilizado o número de
multiplicações e adições aritméticas como uma medida do custo computacional, pois o
número de multiplicações e adições está relacionado diretamente com a velocidade
computacional quando os algoritmos são implementados em computadores digitais de
propósito geral ou em processadores de propósito específico.
DFT
Para
Para o calculo de Fourier pela equação desenvolvida são feitas somas e
multiplicações, sendo o numero de pontos da amostra do sinal, numero fornecido
pelo usuário a fim de se obter uma melhor resolução espectral do sinal.
Dizimação no tempo
Conforme escolhido pelo professor o grupo ficou responsável por desenvolver a
função de dizimação no tempo,
A decomposição e o cálculo da TFD em TFDs sucessivamente menores, explorando-
se tanto a simetria quanto a periodicidade da exponencial complexa
resultam em um aumento drástico na eficiência do cálculo da TFD. Os algoritmos em
que a decomposição é baseada em decompor a sequência em subsequências
sucessivamente menores são chamados de algoritmos de dizimação no tempo.
O princípio da dizimação no tempo é ilustrado de modo conveniente considerando-se
o caso especial de N sendo uma potência inteira de 2, isto é, . Como N é
divisível por dois, podemos considerar a realização do cálculo de separando
em duas sequências de ( ) pontos que consistam em pontos de índice par
e pontos de índice ímpar . Na Figura 1 é mostrada essa
decomposição e também o fato de que a sequência original pode ser recuperada
simplesmente pelo entrelaçamento das duas sequências.
A TFD pode ser representada similarmente em termos das TFDs de e
como se segue no domínio da frequência
Que corresponde à transformada referente à Figura 1, e portanto, pode ser
calculado por
A Equação acima que a TFD de N pontos pode ser computada pelo cálculo das
TFDs de pontos e em em vez de
como normalmente fazemos no caso das TFDs de pontos. A equação
anterior pode ser escrita então como:
O cálculo da Equação acima para todos os valores de exige no máximo
ou multiplicações e adições complexas. É fácil verificar que,
para , o total será menor do que A Equação 7 corresponde a
dividir o cálculo original de pontos em dois cálculos de TFD de pontos. Se
for par, como acontece quando é igual a uma potência de 2, então poderemos
considerar o cálculo de cada uma das TFDs de pontos na Equação 7 quebrando
cada uma das somas nessa equação em duas TFDs de pontos, que então
seriam combinadas para produzir as TFDs de pontos.
A Figura 2 mostra o diagrama de fluxo completo para o cálculo da TFD de 8 pontos,
utilizando as equações e propriedades mencionadas anteriormente.
Para o caso mais geral, mas ainda sendo N uma potência de 2, procederíamos
decompondo as transformadas de pontos em transformadas de pontos e
continuaríamos até ficarmos apenas com transformadas de 2 pontos. Isso requer
estágios de cálculo.
Anteriormente, encontramos que na decomposição original de uma transformada de
pontos em duas transformadas de pontos, o número de multiplicações e adições
complexas exigido era . Quando as transformadas de pontos são
decompostas em transformadas de (N=4) pontos, o fator de é substituído por
de modo que o cálculo total requer então
multiplicações e adições complexas. Se isso pode ser feito no máximo
vezes, de modo que, depois de executar essa decomposição tantas vezes
quanto possível, o número de multiplicações e adições complexas é igual a
.
O cálculo no diagrama de fluxo da Figura 2 pode ser reduzido ainda mais se
explorarmos a simetria e a periodicidade dos coeficientes . O resultado é conhecido
como operação borboleta que requer uma adição complexa e uma subtração
complexa, mas apenas uma multiplicação complexa em vez de duas. A Figura 3
mostra isso. Em particular, o número de multiplicações complexas foi reduzido por um
fator 2 em relação ao número na Figura 2.
Implementação MyDFT
function [Y]=MyDFT(X)
N=length(X);%identifica o tamanho do sinal
Wn = exp(-j*(2*pi/N));%Define o valor de Wn
Y = zeros(1,N);%inicia um vetor de zeros
for k=1:N%inicia o somatorio para cada sinal ponto a ponto
soma=0;%inicia o valor em zero
for n=1:N%somatorio
soma=soma+X(n)*Wn^((k-1)*(n-1));%somatorio dos termos da serie
end
Y(k)=soma;%define o valor no ponto
end
end
Para efetuar a transformada de Fourier é feito uma soma de todos os termos
multiplicados ponto a ponto por um fator , sendo o indicie do termo da
DFT e o índice do termo do sinal.
Implementação MyFFT_Dizimação no tempo
function [ X ] = MyFFT_DecTempo( x )
%Para elementos que nao sao potencias de 2 adiciona-se zeros
n=length(x);
p=ceil(log2(n));
x=[x zeros(1,(2^p)-n)];
%Tamanho do elemento ajustado
n=length(x);
%Constante Wn
Wn=exp((-2*pi*1i)/n);
w=1;
% Caso base da recursao
if n==1
X=x;
else
par = x(1:2:length(x)); % Elementos que compoem os indices
pares do sinal
impar = x(2:2:length(x)); % Elementos que compoem os indices
impares do sinal
Xpar=MyFFT_DecTempo(par);
Yimpar=MyFFT_DecTempo(impar);
for i=1:(n/2)
X(i) = Xpar(i)+Yimpar(i)*w;
X(i+n/2)= Xpar(i) - Yimpar(i)*w;
w=w*Wn; %Armazena valores de Wn para cada operação.
end
end
end
Este código é uma implementação recursiva da dizimação no tempo, para cada
chamada da função, o vetor recebido é dividido em dois, um com índice par e outro
com impares, que são passados recursivamente para a função. O caso base ocorre
quando os vetores contem apenas um elemento. Durante esse processo são
calculadas as funções e .
Para a FFT com dizimação no tempo são realizadas multiplicações
complexas e adições complexas.
Primeiro sinal
Dado um sinal com 2048 pontos, sendo 2048 não 2000 devido a função desenvolvida
trabalhar na potencia de base 2, o sinal é a soma de três senóides com frequências de
15Hz, 45Hz e 70Hz, executando a transformada de Fourier pelos dois métodos
desenvolvidos e pela função FFT do próprio MatLab, obtemos o seguinte resultado.
Na figura é possível observar o sinal e os módulos das transformadas pelos métodos
desenvolvidos e pelo método próprio do software. As transformadas foram colocadas
em figuras diferentes, pois a magnitude atingiu erro aproximadamente zero, a
sobreposição do mesmo não permitiria a visualização de todas elas.
Ao contrario da magnitude a fase houve uma diferença nos resultados, para o caso da
MyDFT, sinal em vermelho, a fase foi aproximadamente zero, então divergiu muito das
demais.
Segundo Sinal
Para a segunda parte foi criado um sinal, que é a soma de duas senóides uma com
frequência de 10Hz e outra com frequência de 25Hz, mas foi calculado para diversos
números de amostras e analisado as comparações da magnitude entre as funções, a
fase entre as funções e o tempo de execução das funções.
As Próximas figuras mostram as comparações para cada numero de amostras.
2 amostras:
4 amostras:
8 amostras:
16 amostras:
32 amostras:
64 amostras:
128 amostras:
256 amostras:
512 amostras:
1024 amostras:
2048 amostras:
Como observado nas figuras, um numero muito baixo de amostras retornam como
resultado valores que não permite a interpretação dos dados do sinal, se tornando
impossível identificar a frequência do sinal e a amplitude a cada frequência no sinal.
Posteriormente foi Verificado o tempo de execução de cada função, sendo em
vermelho a MyDFT, em verde a FFT por dizimação no tempo e em azul a FFT do
MatLab.
Na Tabela abaixo é possível verificar o tempo em segundos.
Comparação Comparação
Amostras MyFFT MyDFT FFT
MyFFtxFFT MyDFTxFFT
2 0,172000 0,618000 0,000000 -1,35E-29 0,00E+00
4 0,127000 0,656000 0,000000 4,26E-28 0,00E+00
8 0,133000 1,352000 0,001700 2,32E-25 3,85E-28
16 0,118000 1,525000 6,641000 -9,75E-25 1,37E-28
32 0,223000 5,370000 0,001600 5,86E-24 8,02E-28
64 0,157000 0,002200 0,003700 -1,61E-21 -3,68E-26
128 0,187000 0,011600 0,006900 6,11E-21 1,30E-25
256 0,278000 0,055400 0,013900 -3,77E-21 -6,96E-25
512 0,438000 0,255400 0,027800 -6,71E-17 -9,90E-23
1024 0,572000 1,151600 0,055400 -5,36E-16 -1,11E-21
2048 0,524000 5,076200 0,114900 -5,38E-14 -7,52E-22
Resultados e Conclusão
Como é possível observar na tabela, o método MyDFT é mais simples mas demora
mais e requer mais custo computacional conforme aumenta o numero de amostras, a
função MyFFT, requer menos tempo na execução mas a FFT é a mais eficiente
comparada as outras funções.
O número de Amostras é importante para uma boa resolução espectral, o uso de um
número muito pequeno de amostras pode atrapalhar muito a amostragem e o
desenvolvimento do experimento.
Anexo
clc;
clear all;
close all;
N=2048;
fs=2048;
t=0:1/N:1-1/N;
vf=0:1:N-1;
X1=sin(2*pi*t*15);
X2=sin(2*pi*t*45);
X3=sin(2*pi*t*70);
X=X1+X2+X3;
Y1=fft(X);
Y2=MyDFT(X);
Y3=MyFFT_DecTempo(X);
F1=unwrap(angle(Y1));
F2=unwrap(angle(Y2));
F3=unwrap(angle(Y3));
figure
subplot(2,2,1);
plot (t,X);
xlabel('Tempo');ylabel('Amplitude');title('Sinal de Entrada');
subplot(2,2,2);
plot (vf,abs(Y1)/N,'b');
xlabel('Frequencia');ylabel('Magnitude');title('Transformada de
Fourier,FFT')
xlim([0 N-1]);
subplot(2,2,3);
plot (abs(Y2)/N,'r');
xlabel('Frequencia');ylabel('Magnitude');title('Transformada de
Fourier,MyDFT')
xlim([0 N-1]);
subplot(2,2,4);
plot (abs(Y3)/N,'g');
xlabel('Frequencia');ylabel('Magnitude');title('Transformada de
Fourier,MyFFT Dizimação no tempo')
xlim([0 N-1]);
figure
subplot(3,1,1);
plot (t,X);
xlabel('Tempo');ylabel('Amplitude');title('Sinal de Entrada');
subplot(3,1,2);
plot(vf,F1);
xlabel('Frequencia');ylabel('Fase');title('Resposta em Fase');
xlim([0 N-1]);
subplot(3,1,3);
hold on;
plot (vf,F1,'b');
plot (F2,'r');
plot (F3,'g');
xlabel('Frequencia');ylabel('Fase');title('Comparação Das Respostas em
Fase');
xlim([0 N-1]);
%% fase e frequencia
tempos=zeros(5,11);
vetor=zeros(1,11);
for i=1:1:11
N = 2^i;
vetor(i)=N;
vf=0:1:N-1;
fs=512;
t=0:1/fs:N/fs-(1/fs);
S1 = sin(2*pi*t*10);
S2 = sin(2*pi*t*25);
S = S1+S2;
tic;
Y1=fft(S);
tempos(1,i)=toc;
tic;
Y2=MyDFT(S);
tempos(2,i)=toc;
tic;
Y3=MyFFT_DecTempo(S);
tempos(3,i)=toc;
tempos(4,i)=erro(Y1,Y2);
tempos(5,i)=erro(Y1,Y3);
figure;
subplot(2,1,1);
plot (vf,abs(Y1)/N,'b');
hold on;
plot (vf,abs(Y2)/N,'r');
hold on;
plot (vf,abs(Y3)/N,'g');
xlim([0 N-1]);
xlabel('Frequencia');ylabel('Magnitude');title('Comparação Das
Respostas em Magnitude');
subplot(2,1,2);
plot (vf,unwrap(angle(Y1)),'b');
hold on;
plot (vf,unwrap(angle(Y2)),'r');
hold on;
plot (vf,unwrap(angle(Y3)),'g');
xlim([0 N-1]);
xlabel('Frequencia');ylabel('Fase');title('Comparação Das Respostas em
Fase');
end;
figure
plot(vetor,tempos(1,:),'b');
hold on;
plot(vetor,tempos(2,:),'r');
hold on;
plot(vetor,tempos(3,:),'g');
xlim([0 2048]);
xlabel('Numero de amostras');ylabel('Tempo');title('Comparação Do
Tempo de Execução');