Slide - Ransformada Rápida de Fourier - Fft - 2024-1

ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 1/19 Processamento Digital de Sinais ELF51 - EL66D - Engenharia Eletrônica Transformada rápida de Fourier - FFT Prof. Daniel R. Pipa danielpipa@utfpr.edu.br ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 2/19 Introdução A DFT (Discrete Fourier Transform) X[k] = N−1 󳕗 n=0 x[n]Wkn N DFT ←−→ x[n] = 1 N N−1 󳕗 k=0 X[k]W−kn N onde WN = e− j 2π N Custo computacional do cálculo direto Tem ordem O(N2) (N2 multiplicações complexas e N(N − 1) adições). Cálculo eficiente É possível reduzir o custo para O(N log2 N). Tais transformadas são chamadas de FFT (Fast Fourier Transform) - Transformada Rápida de Fourier. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 3/19 Custo computacional DFT versus FFT Para valores elevados de N, a vantagem da FFT é significativa. Para N = 1024, a FFT é 100 vezes mais rápida. 0 5 10 15 20 25 30 35 0 200 400 600 800 1000 1200 N O (N ) Computational cost DFT FFT ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 4/19 Analogia com polinômios Seja o poliônomio (ATENÇÃO: variável w e coeficientes xi) X(w) = x0 + x1w + x2w2 + x3w2 + x4w4 + x5w5. Define-se Xe(w) = x0 + x2w + x4w2 Xo(w) = x1 + x3w + x5w2 e X(w) pode ser escrito como X(w) = Xe(w2) + wXo(w2) pois X(w) = Xe(w2) + wXo(w2) = x0 + x2w2 + x4(w2)2 + w 󰀓 x1 + x3w2 + x5(w2)2󰀔 ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 5/19 Analogia com polinômios 2 Para se avaliar o polinômio em pares de pontos wi e −wi X(wi) = Xe(w2 i ) + wiXo(w2 i ) X(−wi) = Xe((−wi)2) − wiXo((−wi)2) Portanto X(wi) = Xe(w2 i ) + wiXo(w2 i ) X(−wi) = Xe(w2 i ) − wiXo(w2 i ) Reaproveitamento de cálculos Para se avaliar o polinômio em −wi, pode-se reaproveitar o cálculo de Xe(w2) e Xo(w2) feito para wi. Simetria Isso é possível pela simetria: Xe((wi)2) = Xe((−wi)2). ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 6/19 Derivação algébrica - Decimação no Tempo Partindo-se da DFT de N pontos, k = 0, 1, · · · , N − 1, X[k] = N−1 󳕗 n=0 x[n]Wkn N = x[0] 󰀓 Wk N 󰀔0 +x[1] 󰀓 Wk N 󰀔1 +x[2] 󰀓 Wk N 󰀔2 +· · ·+x[N −1] 󰀓 Wk N 󰀔 N−1 (ATENÇÃO: polinômio com variável Wk N e coeficientes x[n]). Dividindo o somatório em dois X[k] = N 2 −1 󳕗 n=0 x[2n]Wk(2n) N + N 2 −1 󳕗 n=0 x[2n + 1]Wk(2n+1) N = N 2 −1 󳕗 n=0 xe[n]Wk(2n) N + Wk N N 2 −1 󳕗 n=0 xo[n]Wk(2n) N com xe[n] = x[2n] e xo[n] = x[2n + 1]. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 7/19 Derivação algébrica - Decimação no Tempo cont. Usando a identidade Wk(2n) N = e− j 2π N k(2n) = e − j 2π N 2 kn = Wkn N 2 obtém-se X[k] = N 2 −1 󳕗 n=0 xe[n]Wkn N 2 + Wk N N 2 −1 󳕗 n=0 xo[n]Wkn N 2 = Xe[k] + Wk N Xo[k] com Xe[k] = N 2 −1 󳕗 n=0 xe[n]Wkn N 2 e Xo[k] = N 2 −1 󳕗 n=0 xo[n]Wkn N 2 DFT’s com a metade do tamanho As fórmulas anteriores descrevem duas DFT’s com metade do tamanho da original (Xe[k] e Xo[k] têm tamanho N/2). ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 8/19 Derivação algébrica - Decimação no Tempo cont. Para k = 0, 1, . . ., N/2 − 1, usa-se X[k] = Xe[k] + Wk N Xo[k]. Para k = N/2, . . ., N − 1, nota-se que W k+ N 2 N = Wk NW N 2 N = Wk Ne − j✄2π ✚ N ✚ N ✄2 = −Wk N e explora-se a periodicidade de Xe[k] e Xo[k], já que se tratam de DFT’s de N/2 Xe [k + N/2] = Xe[k] Xo [k + N/2] = Xo[k] Finalmente, X[k] = Xe[k] + Wk N Xo[k], k = 0, 1, · · · , N/2 − 1 X [k + N/2] = Xe[k] − Wk N Xo[k], k = 0, 1, · · · , N/2 − 1 Reaproveita-se Xe[k] e Xo[k] já calculados para X[k], k = 0, . . ., N/2, para X[k + N/2], k = 0, . . ., N/2. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 9/19 Derivação algébrica - Decimação no Tempo cont. Redução do Custo Computacional ◮ Originalmente, N2. ◮ Agora, duas DFT’s (Xe[k] e Xo[k]) mais multiplicações por termo fatorado. 2 󰀕 N 2 󰀖2 + N = N2 + N 2 ◮ Comparando lim N→∞ N2+N 2 N2 = lim N→∞ N2 + N 2N2 = lim N→∞ N + 1 2N = lim N→∞ N 2N = 1 2. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 10/19 Derivação algébrica - Decimação no Tempo cont. Algoritmo recursivo Para se calcular Xe[k] e Xo[k], pode-se aplicar as ideias acima recursivamente. Xe[k] = Xee[k] + Wk N/2Xeo[k], k = 0, 1, · · · , N/4 − 1 Xe [k + N/4] = Xee[k] − Wk N/2Xeo[k], k = 0, 1, · · · , N/4 − 1 Xo[k] = Xoe[k] + Wk N/2Xoo[k], k = 0, 1, · · · , N/4 − 1 Xo [k + N/4] = Xoe[k] − Wk N/2Xoo[k], k = 0, 1, · · · , N/4 − 1 De maneira genérica Xi[k] = Xie[k] + Wk LXio[k], k = 0, 1, · · · , L/2 − 1 Xi [k + L/2] = Xie[k] − Wk LXio[k], k = 0, 1, · · · , L/2 − 1 onde i pode ser e ou o, e L é o tamanho de Xi[k]. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 11/19 Derivação algébrica - Decimação no Tempo cont. até se chegar a uma DFT de 1 ponto. A DFT de 1 ponto é simplesmente X[0] = x[0]W0 1 = x[0] Etapas da FFT recursiva 1. Shuffling (embaralhamento): divide-se sucessivamente a sequência x[n] até N sequências de 1 amostra (DFT’s triviais) 2. Merging (re-combinação): combinam-se sucessivamente as DFT’s até obter-se X[k]. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 12/19 Exemplo Calcular a DFT de x[n] = {7, 5, 9, 3} pelo método de FFT recursiva. Shuffling xe = {7, 9} xo = {5, 3} xee = {7} xeo = {9} xoe = {5} xoo = {3} DFT’s triviais Xee = {7} Xeo = {9} Xoe = {5} Xoo = {3} Merging ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 13/19 Exemplo cont. L = 2 Xe[k] = Xee[k] + Wk 2 Xeo[k], k = 0 Xe [k + 1] = Xee[k] − Wk 2 Xeo[k], k = 0 Xo[k] = Xoe[k] + Wk 2 Xoo[k], k = 0 Xo [k + 1] = Xoe[k] − Wk 2 Xoo[k], k = 0 W0 2 = 1 Xe[0] = Xee[0] + W0 2 Xeo[0] = 7 + 1 × 9 = 16 Xe [1] = Xee[0] − W0 2 Xeo[0] = 7 − 1 × 9 = −2 Xo[0] = Xoe[0] + W0 2 Xoo[0] = 5 + 1 × 3 = 8 Xo [1] = Xoe[0] − W0 2 Xoo[0] = 5 − 1 × 3 = 2 ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 14/19 Exemplo cont. L = 4 = N X[k] = Xe[k] + Wk 4 Xo[k], k = 0, 1 X [k + 2] = Xe[k] − Wk 4 Xo[k], k = 0, 1 W0 4 = 1, W1 4 = e−j 2π 4 = −j X[0] = Xe[0] + W0 4 Xo[0] = 16 + 1 × 8 = 24 X[1] = Xe[1] + W1 4 Xo[1] = −2 − j × 2 = −2 − j2 X[2] = Xe[0] − W0 4 Xo[0] = 16 − 1 × 8 = 8 X[3] = Xe[1] − W1 4 Xo[1] = −2 + j × 2 = −2 + j2 ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 15/19 Custo computacional da FFT FFT recursiva X[k] = Xe[k] + Wk N Xo[k], k = 0, 1, · · · , N/2 − 1 X [k + N/2] = Xe[k] − Wk N Xo[k], k = 0, 1, · · · , N/2 − 1 ◮ Cada etapa de merging dobra o tamanho do resultado. Em r etapas, são geradas sequências de 2r amostras. ◮ Para o resultado final, é preciso N = 2r amostras ou r = log2 N etapas de merging. ◮ Cada etapa de merging requer N multiplicações complexas. ◮ O custo final é O(N log2 N). Exemplo Python ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 16/19 Algoritmo de Goertzel Motivação A FFT, embora seja rápida, calcula a DFT para todos os valores de k. Algumas aplicações requerem o cálculo da DFT somente para algumas frequências (p. ex. detecção de senoides puras). Algoritmo de Goertzel O algoritmo de Goertzel é mais rápido que a FFT se forem necessárias apenas M < log2 N pontos da DFT. Sabendo que W−kN N = e j 2π N kN = 1 tem-se X[k] = W−kN N N−1 󳕗 n=0 x[n]Wkn N = N−1 󳕗 n=0 x[n]W−k(N−n) N ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 17/19 Algoritmo de Goertzel cont. Por exemplo, para N = 4 X[k] = 3 󳕗 n=0 x[n]W−k(4−n) 4 = x[3]W−k 4 + x[2]W−2k 4 + x[1]W−3k 4 + x[0]W−4k 4 = W−k 4 󰀋 x[3] + W−k 4 󰀋 x[2] + W−k 4 󰀋 x[1] + W−k 4 x[0] 󰀌󰀌󰀌 Generalizando e fazendo yk[−1] = 0, tem-se a definição recursiva yk[n] = x[n] + W−k N yk[n − 1], 0 ≤ n ≤ N X[k] = yk[N] Filtro IIR A recursividade acima parece um filtro IIR! ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 18/19 Algoritmo de Goertzel cont. A função de transferência é Hk(z) = 1 1 − W−k N z−1 = 1 − Wk N z−1 󰀃1 − W−k N z−1󰀄 󰀃1 − Wk N z−1󰀄 = 1 − Wk N z−1 1 − 2 cos(2πk/N)z−1 + z−2 = 󰀗 1 1 − 2 cos(2πk/N)z−1 + z−2 󰀘 󰀅 1 − Wk N z−1󰀆 o que resulta em vk[n] = 2 cos(2πk/N)vk[n − 1] − vk[n − 2] + x[n], 0 ≤ n ≤ N X[k] = vk[N] − Wk Nvk[N − 1] com condições iniciais vk[−1] = vk[−2] = 0. Exemplo Python. ELF51 - EL66D - PDS Prof. Daniel R. Pipa Transformada rápida de Fourier - FFT Introdução Decimação no Tempo Exemplos Algoritmo de Goertzel 19/19 Custos computacionais Multiplic. reais Multiplic. complexas DFT 0 O(N2) Goertzel 1 0 O(MN) Goertzel 2 O(MN) O(M) FFT 0 O(N log2 N) Outras vantagens ◮ Memória: Por causa de recursividade, não é necessário armazenar N valores complexos. ◮ Eficiência: Em aplicações de tempo real, os cálculos podem começar assim que a primeira amostra x[0] estiver disponível, dispensando encher o “buffer”.