Szybkie metody obliczania DFT i ich zastosowania
Autorzy: Krzysztof Duda, Tomasz P. Zieliński
Pliki: rozdzial_05.zip
Przykład 5.1
W przykładzie 5.1 (plik fft_dit.m) jest zamieszczona implementacja w języku MATLAB algorytmu FFT z podziałem w czasie.% Przykład 5.1 - Algorytm FFT z podziałem w czasie (DIT) function Xw=fft_dit(xt); N=length(xt); %% uzupełnienie zerami do długości 2^v v=ceil(log2(N)); %liczba bitów xt=[xt zeros(1,2^v-N)]; N=length(xt); %% odwrócenie kolejności bitów indeksowanie od 0 xo= zeros(size(xt)); %inicjalizacja zmiennej b = 2.^(v-1:-1:0); %wagi binarne for k=0:N-1; ind=k; ko=zeros(1,v); for k1=0:v-1 if ind-b(k1+1)>=0 ko(k1+1)=1; ind=ind-b(k1+1); end end ind_o=sum(fliplr(ko).*b); xo(ind_o+1)=xt(k+1); end %% algorytm FFT X=zeros(2,N); %pamieć dla danych i wyników X(1,:)=xo; %dane WN=exp(-j*2*pi/N); for k=0:v-1 %petla po etapach M=2^k; %liczba motylków w bloku for k1=0:N/2/M-1; %pętla po blokach for k2=0:M-1; %pętla wewnątrz bloku %ustalenie indeksów p=k1*N/2^(v-k-1)+k2; q=p+M; r=2^(v-k-1)*k2; %obliczenia motylkowe X(2,p+1)=X(1,p+1)+WN^r*X(1,q+1); X(2,q+1)=X(1,p+1)-WN^r*X(1,q+1); end end X(1,:)=X(2,:); %obliczenia z podstawianiem end Xw=X(1,:);Przykład (plik przyklad_5_1.m) ilustrujący użycie funkcji fft_dit()
clear all, close all, clc
N=32; %liczba próbek sygnału
n = 0:N-1;
x=cos(5.5*(2*pi/N)*n)+0.1*randn(1,N); %przykładowy sygnał testowy
X1=fft(x); k1=0:N-1;
X2=fft_dit(x); k2=0:length(X2)-1;
if length(X2)==N
blad = max(abs(X1-X2))
end
figure
plot(n, x,'.-b')
xlabel('n')
ylabel('x_n')
title('Sygnal testowy')
box on, axis tight
figure, hold on
subplot(3,1,1), hold on
plot(k1, abs(X1),'.-b')
plot(k2, abs(X2),'o-r')
xlabel('k')
title('Modul')
legend('fft','fft_dit')
box on, axis tight
subplot(3,1,2), hold on
plot(k1, real(X1),'.-b')
plot(k2, real(X2),'o-r')
xlabel('k')
title('Czesc rzeczywista')
legend('fft','fft_dit')
box on, axis tight
subplot(3,1,3), hold on
plot(k1, imag(X1),'.-b')
plot(k2, imag(X2),'o-b')
xlabel('k')
title('Czesc urojona')
legend('fft','fft_dit')
box on, axis tight
Przykład 5.2
Wprzykładzie 5.2 (plik przyklad_5_2.m) przedstawiono implementacja w języku MATLAB metody obliczania widma DFT sygnału rzeczywistego o długości 2N, z użyciem jednego N punktowego DFT.close all, clear all, clc N = 4; %liczba probek x = rand(1,2*N); %sygnał testowy o długości 2N WN=exp(-j*2*pi*(0:(2*N)/2-1)/(2*N)); x1=x(1:2:end); % N próbek o indeksach parzystych x2=x(2:2:end); % N próbek o indeksach nieparzystych x12=x1+j*x2; % 'sztuczny' sygnal zespolony X12=fft(x12); % FFT N-punktowe X0 =X12(1); X12(1)=[]; X1 = [real(X0) (X12+fliplr(conj(X12)))/2 ]; %DFT próbek o indeksach parzystych X2 = [imag(X0) (X12-fliplr(conj(X12)))/(2*j)]; %DFT próbek o indeksach nieparzystych X = [X1+WN.*X2]; %podział w częstotliwości X = [X X1(1)-X2(1) conj(X(end:-1:2))]; %% porównanie z Matlabem Xf=fft(x); [X(:)-Xf(:)]
Przykład 5.3
W przykładzie 5.3 (plik przyklad_5_3.m) przedstawiono implementację algorytmu mSDFT w języku MATLAB.close all, clear all, clc
x = rand(1,100); %sygnał testowy
M = 30; %długość DFT
k = 3; %prążek DFT, który ma być wyznaczany
W = exp(-j*2*pi/M); % twiddle factor
%% Wartosci początkowe
xmr = 1;
X0M = 0;
xn_M= 0;
fifo= zeros(1,M);
for n=1:length(x);
%% bufor FIFO
fifo= [fifo(2:M) x(n)];
%% mSDFT
X1M = X0M+xmr*( x(n)-xn_M );
X0M = X1M;
if mod(n,M)
xmr=xmr*W^(k); %ciąg modulujący
X1M = X1M*conj(xmr);%korekcja fazy
else
xmr=1;
end
%% nowe xn_M
xn_M = fifo(1);
end
%% błąd amplitudy i fazy w odniesieniu do blokowego DFT
Xk = fft(fifo); Xk=Xk(k+1); %blokowe DFT
Aerr = [abs(X1M)-abs(Xk)] %bład amplitudy
Perr = [angle(X1M)-angle(Xk)] %bład fazy
Przykład 5.4
W przykładzie 5.4 (plik przyklad_5_4.m) przedstawiono implementację algorytmu szybkiego splotu dwóch sygnałów z wykorzystaniem FFT.% Przykład 5.4
close all, clear all, clc
M = 512; % długość filtra
N = 10*M; % długość sygnału
x=randn(1,N); % sygnał 1
h=0.99.^(0:M-1); % sygnał 2, filtr
% Splot liniowy - funkcja Matlaba
y1 = conv(x,h);
% Szybki splot kołowy całego sygnału za pomocą FFT
hz = [ h zeros(1,N-M)]; % uzupełniamy zerami drugi sygnał
% do długości pierwszego, dłuższego
y2 = ifft( fft(x) .* fft( hz ) ); % szybki splot kołowy
blad2 = max(abs(y1(M:N) - y2(M:N))) % błąd, końcowe N-M+1 próbki są poprawne
% Szybki splot liniowy całego sygnału za pomocą FFT
xz = [ x zeros(1,M-1)]; % uzupełniamy dłuższy sygnał M-1 zerami
hzz = [ h zeros(1,N-M) zeros(1,M-1)]; % uzupełniamy krótszy sygnał N-1 zerami
y3 = ifft( fft(xz) .* fft( hzz ) ); % szybki splot liniowy
blad3 = max(abs( y1- y3)) % błąd, wszystkie próbki są poprawne
Przykład 5.5
W przykładzie 5.5 (plik przyklad_5_5.m), kontynuacji przykładu 5.4, przedstawiono implementację różnych wersji algorytmu sekcjonalnego szybkiego splotu overlap-add.% Przykład 5.5
close all, clear all, clc
Nh = 512; Nx = 10*Nh; % długość filtra, długość sygnału
Nb = Nh; % długość bloku próbek sygnału
x=randn(1,Nx); h=0.99.^(0:Nh-1); % sygnał 1 = losowy, sygnał 2 = odp. impulsowa filtra
y1 = conv(x,h); % wynik splotu funkcją Matlaba jako odniesienie
% Metoda czwarta: sekcjonowany szybki splot overlap-add #1 (dodawanie sygnałów na wyjściu)
Lb = ceil(Nx/Nb); % liczba bloków
hzz = [h zeros(1,Nb-Nh) zeros(1,Nh)]; % wypelnianie h zerami
Hz = fft(hzz); % FFT odpowiedzi impulsowej
y4 = zeros(1,Nh-1); % przygotowanie wektora y4
for k = 1:Lb
bx = x((k-1)*Nb+1:k*Nb); % wycinanie bloku Nb próbek sygnału x[n]
bxz = [bx zeros(1,Nh)]; % uzupełnianie zerami do długości Nb+Nh-1
Bxz = fft(bxz); % FFT fragmentu sygnału
by = real( ifft( Bxz.*Hz ) ); % iloczyn widm i IFFT
y4(end-(Nh-2):end) = y4(end-(Nh-2):end)+ by(1:Nh-1); % dodawanie koniec + początek
y4 = [ y4 by(Nh:end-1)]; % wstawianie próbek na koniec wynikowego wektora
end;
blad4 = max(abs(y1-y4)) % błąd
% Metoda piąta: sekcjonowany szybki splot overlap-add #2 (dodawanie widm na wejściu)
Lb = ceil(Nx/Nb); % liczba bloków próbek
hzz = [h zeros(1,Nb-Nh) zeros(1,Nb)]; % wypełnianie h zerami Nb+Nb
Hz = fft(hzz); % FFT odpowiedzi impulsowej
J = (-1).^(0:2*Nb-1); % J = macierz korekcji opóźnienia w widmie X2
X2 = zeros(1,2*Nb); y5 = []; % przygotowanie wektora y5
for k = 1:Lb
x1 = x((k-1)*Nb+1:k*Nb); % wycinanie bloku Nb próbek sygnału x[n]
xx = [x1 zeros(1,Nb)]; % uzupełnianie zerami do długości Nb+Nh
X1 = fft(xx); % FFT
X = X1 + J.*X2; X2 = X1; % dodawanie widm X1 i X2 (z korekcją opóźnienia)
yy = real( ifft( X.*Hz ) ); % iloczyn widm, następnie IFFT
y5 = [ y5 yy(1:Nb)]; % wstawianie próbek na koniec wynikowego wektora
end;
blad5 = max(abs(y1(1:length(y5))-y5)) % błąd
% Metoda szósta: sekcjonowany szybki splot overlap-add #3 z dodatkowym podziałem h
% na bloki (fragmenty) o równej długości
Nb = 32; % arbitralnie przyjęta długość bloku próbek x[n] i h[n]
Lb = ceil(Nx/Nb); % liczba bloków sygnału
Lh = ceil(Nh/Nb); % liczba bloków odpowiedzi impulsowej filtra
for k = 1:Lh % obliczenie macierzy H (widma fragmentów odpowiedzi impulsowej)
H(k,1:2*Nb) = fft( [h((k-1)*Nb+1:k*Nb) zeros(1,Nb)] ); % uzupełnianie zerami, FFT
end
bx = zeros(1,2*Nb); X = zeros(Lh,2*Nb); y6 = [];
for k = 0 : Lb-1
bx = [ bx(Nb+1:2*Nb) x(k*Nb+1:k*Nb+Nb) ]; % pobierz kolejny blok próbek x[n]
X = [ fft(bx); X(1:Lh-1,:) ]; % oblicz jego FFT, zmodyfikuj macierz X
% (w wierszach widma fragmentów sygnału x)
XH = X.*H; % iloczyn elementów macierzy widm X i H
by = real( ifft( sum(XH) ) ); % sumowanie, potem IFFT
y6 = [ y6 by(Nb+1:2*Nb) ]; % wstawianie próbek na koniec wynikowego wektora
end
blad6 = max(abs(y1(1:length(y6))-y6)) % błąd
% Metoda siódma: sekcjonowany szybki splot overlap-add #4 z dodatkowym podziałem h
% na bloki (fragmenty) o nierównej długości - małe opóźnienie sygnału na wyjściu
N = 32; % rozmiar pierwszego, najkrótszego fragmentu h
max_index = Nh/(4*N); % długość najdłuższego fragmentu h (wyrażona w wielokrotności 4N)
H = zeros((log2(max_index)+1)*2, Nh/2); % inicjalizacja macierzy H widm FFT fragmentów h
Lb = size(H,1); % liczba fragmentów (bloków) h
y7 = zeros(1,Nx+2*Nh-1);
% Podział odpowiedzi impulsowej filtra
i=1; k=1;
while i<=Lb
H(k ,1:i*N) = h(2*N*i+1 : 2*N*i+N*i); % bloki wag filtra są wpisywane do macierzy
H(k+1,1:i*N) = h(3*N*i+1 : 3*N*i+N*i); % o wymiarze: Nh/2 x log2(max_index)
i=i*2; k=k+2;
end
% FFT fragmentów odp. impulsowej filtra (w jednym obiegu pętli liczymy dwa widma)
i=1;
for k=1:2:size(H,1)
H(k, 1:2*i*N) = fft(H(k, 1:2*i*N)); % kazdy odcinek jest wypelniony zerami do
H(k+1,1:2*i*N) = fft(H(k+1, 1:2*i*N)); % dl 2*dlugość odcinka
i=i*2;
end;
% Główna pętla
x = [ x zeros(1,Nh-1) ]; % uzupełnienie sygnału zerami
for n = 1 : Nx+Nh-1 % powtarzamy do końca sygnału
i = 1;
% Pierwszy fragment odpowiedzi impulsowej filtra - metoda bezposrednia
if n<2*N
splot = h(n:-1:1)*(x(1:n))'; % splot w dziedzinie czasu
else
splot = h(2*N:-1:1)*(x(n-2*N+1:n))'; % splot w dziedzinie czasu
end
y7(n) = y7(n)+splot;
% Kolejne fragmenty odpowiedzi impulsowej filtra - metoda blokowa
k = 1;
while(i <= max_index)
if mod(n,i*N)==0 % czy mamy odp. liczbę probek (32, 64, 128, ...)
bx=[x(n-i*N+1:n),zeros(1,i*N)]; % wycinamy fragment sygnału
BX = fft(bx); % FFT
XH1 = BX.*H(k, 1:2*(i*N)); % w jednym przebiegu pętli
XH2 = BX.*H(k+1, 1:2*(i*N)); % liczymy dwa sploty
i1 = n+i*N; % gdzie wstawic wynik?
i2 = n+2*i*N; % to samo dla drugiego odcinka
Ny = 2*i*N; % długość odcinka sygnalu
y7(i1+1:i1+Ny) = y7(i1+1:i1+Ny) + real(ifft(XH1));
y7(i2+1:i2+Ny) = y7(i2+1:i2+Ny) + real(ifft(XH2));
k=k+2;
else
break % jeżeli liczba próbek nie dzieli się przez N,
% to nie dzieli się też przez 2*N, 4*N itd.
end
i=i*2;
end
end
y7=y7(1:length(y7)-Nh);
blad7 = max(abs(y1-y7))
Przykład 5.6
W przykładzie 5.6 (plik przyklad_5_6.m) przedstawiono implementację szybkiego algorytmu wyznaczania korelacji wzajemnej (lub własnej kiedy y[n] = x[n]).% Przykład 5.6
clear all; close all;
% Inicjowanie parametrów
N = 256; % liczba próbek obu sygnałów
M = 64; % liczba przesunięć w jedną stronę
L = 1000; % liczba prążków funkcji wzajemnej gęstości widmowej mocy
m = N-M: N+M; % indeksy współczynników liczonych "szybko"
% Generowanie (synteza) sygnałów
fx = 100; fpr = 1000; dt=1/fpr; t=0:dt:(N-1)*dt; tR=-(N-1)*dt:dt:(N-1)*dt;
x1 = rand(1,N); x2 = randn(1,N);
% x1 = sin(2*pi*fx*t); x2=x1;
subplot(211); plot(x1); grid; axis tight; title('Sygnały analizowane');
subplot(212); plot(t,x2); grid; axis tight; xlabel('czas [sek]'); pause
% Obliczenie funkcji korelacji wzajemnej sygnałów procedurą Matlaba
% R1 = xcorr(x1,x2,'biased'); % tak
R1 = conv(x1, conj(x2(end:-1:1)))/N; % lub tak
% Obliczenia "ręczne" funkcji korelacji wzajemnej bezpośrednio z definicji
for k=0:M
R21(k+1)= sum( x1(1+k:N).*x2(1:N-k) ) / N;
R22(k+1)= sum( x1(1:N-k).*x2(1+k:N) ) / N;
end
R2 = [ R22(M+1:-1:2) R21(1:M+1) ];
subplot(211); plot(tR(m),R1(m)); grid; axis tight; title('Funkcja korelacji - Matlab');
subplot(212); plot(tR(m),R2); grid; axis tight; title('Funkcja korelacji - "ręcznie"'); pause
max(abs( R1(m)-R2)), pause
% Obliczenia "szybkie" 2*M+1 środkowych współczynników funkcji korelacji wzajemnej
% z FFT
x1 = [ x1 zeros(1,M) ]; % dodaj M zer na końcu sygnału
x2 = [ x2(M+1:end) zeros(1,M) x2(1:M) ]; % dodaj M zer "wewnątrz" sygnału
X1 = fft(x1); % oblicz FFT (N+M)-punktowe
X2 = fft(x2); % oblicz FFT (N+M)-punktowe
X = X1.*conj(X2); % oblicz X1(k)*conj(X2(k))
R = ifft(X); % oblicz odwrotne FFT (N+M)-punktowe
R3 = real(R(1:2*M+1))/N; % pobierz poprawne wartości, przeskaluj
subplot(111); plot(tR(m),R1(m),'r',tR(m),R3,'b'); grid; axis tight;
title('Metoda "szybka" - Matlab (red), My (blue)'); xlabel('czas [sek]'); pause
max(abs( R1(m)-R3)), pause
% Funkcja wzajemnej gęstości widmowej mocy metodą Blackmana-Tukeya
w = hanning(2*M+1); w=w'; % wybierz okno widmowe, np. Hanninga
Rw = R3 .* w; % wymnóż funkcję autokorelacji z oknem
s = [ Rw(M:2*M+1) zeros(1,L-2*M+1) Rw(1:M+1)]; % symetria dane wejściowe dla FFT
S = abs(fft(s)); % transformacja Fouriera
df = 1/(L*dt); f=0:df:(L-1)*df; k=1:L/2+1;
subplot(111); plot(f(k),S(k),'b'); grid; title('Funkcja gęstości widmowej mocy');
xlabel('f [Hz]'); pause
Przykład 5.7
W przykładzie 5.7 (plik przyklad_5_7.m) przedstawiono implementację szybkiego obliczania DTFT (Discrete-Time Fourier Transform) z wykorzystaniem FFT.% Przykład 5.7
clear all; close all;
% Inicjowanie wartości parametrów
N = 256; % liczba próbek sygnału
M = 32; % liczba prążków widma
fpr = 128; % częstotliwość próbkowania
fd = 10; fg = 20; % częstotliwości dolna i górna
df=(fg-fd)/(M-1); f=fd:df:fg; % interesujące nas częstotliwości
x = rand(1, N); % sygnał analizowany, szum z przedziału [0, 1]
NM1 = N+M-1;
A = exp( -j*2*pi * fd/fpr ); % punkt startowy
W = exp( -j*2*pi * ((fg-fd)/(2*(M-1))/fpr) ); % krok
y1 = zeros(1,NM1); k=0:N-1; y1(k+1)= ((A*W.^k).^k) .* x(k+1); % inicjowanie y1
k = [ 0 : M-1, -(N-1):1:-1]; y2 = W.^(-k.^2); % inicjowanie y2
%k = [ 0:-1:-(N-1), M-1:-1:1]; k = [ k(1) fliplr(k(2:end))]; y2 = W.^(-k.^2);
% inicjowanie y2
Y1 = fft(y1); % # algorytm szybkiego splotu kołowego
Y2 = fft(y2); % # sygnałów y1 i y2
Y = Y1.*Y2; % #
y = ifft(Y)/(N/2); % #
k=0:M-1; XcztN(k+1) = y(k+1) .* (W.^(k.^2)); % korekcja fazowa wyniku splotu
n = 0:N-1; Xref=x*exp(-j*2*pi*n(:)*f/fpr); % obliczenia na podstawie definicji DTFT (2.14)
Xref = Xref/(N/2); % skalowanie
blad = max(abs( XcztN - Xref )), % błąd w stosunku do obliczeń z definicji
subplot(2,1,1);
plot(f,real(XcztN),'r.-'); grid on; hold on;
plot(f,real(Xref),'bo-');
subplot(2,1,2);
plot(f,imag(XcztN),'r.-'); grid on; hold on;
plot(f,imag(Xref),'bo-'); pause