Zapraszam do lektury drugiego wpisu, który dotyczy próbkowania sygnałów oraz dyskretnego przekształcenia Fouriera DFT w MATLABie. Tym razem przyjrzymy się konsekwencjom próbkowania oraz cechom widma częstotliwościowego. Po lekturze tego wpisu pojęcie symetrii DFT czy przecieku nie będzie dla Was tajemnicą!

Z poprzedniego postu wiemy, że w świecie sygnałów dyskretnych występuje niejednoznaczność związana z próbkowaniem sygnałów ciągłych: ten sam zbiór próbek może reprezentować nieskończoną liczbę sygnałów o częstotliwościach: f0 + k*fpr. Ponadto rozważania zawarte w poprzednim wpisie skierowały naszą uwagę na ważny problem: z jaką częstotliwością należy próbkować sygnał ciągły, aby otrzymany zbiór próbek zachował jego wartość informacyjną? Oraz jakie są dalsze konsekwencje próbkowania sygnałów?

DFT - Dyskretna transformacja Fouriera

Sygnał ciągły x(t) ma swoją ciągłą transformatę Fouriera. Kiedy próbkujemy taki sygnał to otrzymujemy ciąg próbek x(n) o liczności n=1,2,…N. Jeżeli dla sygnału spróbkowanego obliczymy DFT, to otrzymamy spróbkowaną aproksymacje (przybliżenie) ciągłej transformaty Fouriera sygnału ciągłego. Im więcej próbek posiada analizowany fragment sygnału (czyli im większa częstotliwość próbkowania), tym lepiej DFT przybliża prawdziwe widmo sygnału ciągłego. Weźmy następujący przykład. Sygnał ciągły zawiera dwie harmoniczne o częstotliwościach 5Hz i 15.5Hz i amplitudach odpowiednio 1 i 0.5. Sygnał próbkujemy z częstotliwością fpr=50Hz i zbieramy N=50 próbek sygnału. Do wyznaczenia widma sygnału używamy algorytmu FFT, według poniższego schematu (cały program jest na końcu tego postu).

%% Generacja sygnałów dyskretnych (próbkowanie)
fpr=50; dt=1/fpr;       % Częstotliwość i okres próbkowania
t0=1; N=t0/dt;          % Czas obserwacji sygnałów i liczba próbek
t=(0:N-1)*dt;           % Wektor czasu
df=fpr/N; F=df*(0:N-1); % Wektor częstotliwości do FFT

…..	% tu generujemy sygnały

%% Obliczanie DFT z wykorzystaniem algorytmu FFT i rysunki
X = abs(fft(x)); 
 
figure;
stem(F, X); 

Zobaczmy na rysunki:

1. Sygnał ciągły narysowany w funkcji czasu oraz jego próbki:

2. Widmo sygnału ciągłego (prawdziwe – tego nie znamy, ale chcemy policzyć przez FFT) :

3. Widmo sygnału spróbkowanego (otrzymane na drodze obliczeń z algorytmem FFT):

Już na pierwszy rzut oka widać, że po obliczeniu DFT za pomocą algorytmu FFT coś jest nie tak. Obraz z ostatniego rysunku nie odpowiada rysunkowi transformaty Fouriera sygnału ciągłego. Po pierwsze prążek przy częstotliwości 15.5Hz wydaje się być rozmyty. Po drugie w widmie występuje więcej częstotliwości niż jest ich naprawdę, a widmo wydaje się być symetryczne względem połowy częstotliwości próbkowania. Po trzecie wartości harmonicznych nie odpowiadają wartościom 1 i 0.5. Gdyby więc wprost interpretować ostatni rysunek, to stwierdzilibyśmy, że w sygnale występują częstotliwości 5Hz, 45Hz o wartościach 25, oraz szereg częstotliwości bliskich 15Hz i 35Hz.. Nic bardziej mylnego! Co się tu zatem stało?

Powyższe obserwacje są konsekwencją próbkowania sygnału oraz cech dyskretnej transformaty Fouriera. Przyjrzyjmy się temu z bliska.

Rozdzielczość DFT

Mając N próbek sygnału w dziedzinie czasu, wyznaczamy widmo sygnału również w N punktach, równomiernie rozłożonych na osi częstotliwości od 0 do fpr. Liczba próbek N jest więc ważnym parametrem gdyż określa jaka jest rozdzielczość DFT w dziedzinie częstotliwości. Częstotliwości dla których jest wyznaczana DFT są opisane wyrażeniem: f_DFT = m*fpr/N, gdzie m = 0, 1, 2, …N.

Symetria DFT

Z rysunku widać, że widmo sygnału jest symetryczne  względem połowy częstotliwości próbkowania. Faktycznie, z teorii wiemy, że dla rzeczywistych sygnałów dyskretnych x(n), wartości DFT dla częstotliwości m>(N/2), są „nadmiarowe”, gdyż prążek dla m-tej częstotliwości będzie miał taką samą wartość jak dla częstotliwości (N-m). Innymi słowy, jeżeli dla sygnału x(n) wyznaczymy N punktową DFT, to otrzymamy N prążków reprezentujących występujące w sygnale częstotliwości, ale tylko pierwsze N/2 jest niezależnych. Dlatego licząc i rysując widmo sygnału, np. za pomocą algorytmu FFT w MATLABie, robimy to dla pierwszych N/2 wartości częstotliwości, czyli 0<= m <= (N/2)-1. Wartości widma powyżej N/2 nie wnoszą żadnej dodatkowej informacji, a są jedynie „odbiciem” prążków w paśmie podstawowym.

Przeciek widma

Wiemy już, że DFT jest wyznaczana dla dyskretnych wartości częstotliwości. Bardzo ważną konsekwencją tego, że DFT „próbkuje” widmo sygnału ciągłego jest to, że jest ono poprawnie wyznaczane tylko dla harmonicznych sygnału x(n) o częstotliwościach odpowiadających częstotliwościom w których jest obliczane DFT, czyli dla całkowitych wielokrotności  fpr/N. Jeżeli sygnał zawiera częstotliwość np. 1.5*fpr/N, czyli pomiędzy częstotliwościami dla których jest liczone DFT m*fpr/N, to energia związana z tą częstotliwością pojawi się we wszystkich innych N częstotliwościach obliczonego widma. Zjawisko to nazywa się przeciekiem widma, bo faktycznie.. energia harmonicznej 1.5*fpr/N  jakby „wycieka” na sąsiednie prążki widma. Z taką sytuacją mamy właśnie do czynienia w naszym przykładzie w okolicy częstotliwości 15.5Hz. Jak można się domyśleć nie jest to najlepsza wiadomość, gdyż w większości przypadków w rzeczywistości mamy do czynienia z sygnałami o składowych „pomiędzy” punktami wyznaczania widma. Nasz wynik będzie więc bardziej lub mniej „rozmyty” przez efekt przecieku. Zjawisko to można minimalizować, między innymi, poprzez dobór odpowiednich okien czasowych ale jest temat na osobny artykuł.

Wartości widma DFT

Z przedstawionego przykładu widać, że wartości prążków dla częstotliwości 5Hz i 15.5Hz, są większe niż amplitudy harmonicznych tworzących sygnał x(n). Dlatego należy brać pod uwagę czynnik skalujący, który dla sygnałów rzeczywistych wynosi N/2.

Biorąc pod uwagę powyższe informacje należy zmodyfikować program obliczający DFT w MATLABie w następujący sposób:

%% Obliczanie DFT z wykorzystaniem algorytmu FFT i rysunki
X = (2/N)*abs(fft(x)); 
 
figure;
stem(F(1:N/2), X(1:N/2)); 

Po wykonaniu programu otrzymamy:

Jest więc dużo lepiej! Nie obserwujemy nadmiarowego pasma powyżej połowy częstotliwości próbkowania, a wartości harmonicznych odpowiadają (przynajmniej tej 5Hz) wartościom prawdziwym. Problemem wciąż jest przeciek widma wokół częstotliwości 15.5Hz. Ale to, jak już wspominałem, jest temat na osobny wpis.

Powoli zbliżamy się do odpowiedzi na fundamentalne pytanie o częstotliwość próbkowania sygnałów i jej wpływ na wyniki obliczeń. Tym zajmę się w kolejnej, trzeciej części poświęconej DFT.

Załącznik (program generujący powyższe rysunki):

clear all; close all; clc

%% Generacja sygnałów "ciągłych"  
% Musimy je "symulować" w MATLABie, gdyż w komputerze każdy sygnał jest
% dyskretny. Stąd słowo "ciągły" jest w cudzysłowie - efekt ciągłości 
% uzyskujemy przez ustawienie dużej częstotliwości próbkowania. 
% Robimy to tylko po to, aby móc narysować sygnały "ciągłe" na tle ich
% spróbkowanych odpowiedników. Sygnały "ciągłe" nie biorą udziału w analizie
% częstotliwościowej. 

fpr_cont=10e3; dt_cont=1/fpr_cont;      % wysoka częstotliwość próbkowania
to_cont=1; N_cont=to_cont/dt_cont;      % czas obserwacji sygnałów
t_cont=(0:N_cont-1)*dt_cont;            % wektor czasu "ciągłego"

A1_cont = 1;                    % Amplitudy harmonicznych
A2_cont = 0.5;

f1_cont=5; fi1_cont=0;          % częstotliwości i fazy harmonicznych "ciągłych"
f2_cont=15.5; fi2_cont=0; 

x1_cont = A1_cont*sin(2*pi*f1_cont*t_cont+fi1_cont);    % generacja sygnałów "ciągłych"
x2_cont = A2_cont*sin(2*pi*f2_cont*t_cont+fi2_cont);
x_cont = x1_cont + x2_cont;

%% Generacja sygnałów dyskretnych (próbkowanie)
fpr=50; dt=1/fpr;       % Częstotliwość i okres próbkowania
t0=1; N=t0/dt;          % Czas obserwacji sygnałów i liczba próbek
t=(0:N-1)*dt;           % Wektor czasu
df=fpr/N; F=df*(0:N-1); % Wektor częstotliwości do FFT

% Parametry sygnałów dyskretnych oczywiście takie jak ciągłych
A1 = A1_cont;
A2 = A2_cont;

f1=f1_cont; fi1=fi1_cont;
f2=f2_cont; fi2 = fi2_cont;

% Generacja sygnałów dyskretnych
x1 = A1*sin(2*pi*f1*t+fi1);
x2 = A2*sin(2*pi*f2*t+fi2);
x = x1 + x2;

%% Rysunki
figure; 
plot(t_cont, x_cont)
hold on; plot(t, x, 'ob'); grid on;
title(['Sygnał o składowych: ' num2str(f1) ' oraz ' num2str(f2) ' [Hz]']); 
xlabel('Czas [s]'); ylabel('Amplituda [-]')
set(findobj(gcf,'type','axes'),'FontName','Halvetica','FontSize',12);

%% Prawdziwe widmo
figure; 
stem([f1_cont, f2_cont], [A1_cont, A2_cont]);
title(['Widmo sygnału ciągłego o składowych ' num2str(f1_cont) 'Hz oraz ' num2str(f2_cont) 'Hz'])
xlabel('Częstotliwość [Hz]'); ylabel('Amplituda [-]');
axis([0 50 0 1.5]); grid on;
set(findobj(gcf,'type','axes'),'FontName','Halvetica','FontSize',12);

%% Obliczanie DFT z wykorzystaniem algorytmu FFT i rysunki
X = (2/N)*abs(fft(x)); 

figure;
stem(F(1:N/2), X(1:N/2)); 

grid on;
title(['Widmo sygnału dyskretnego o składowych ' num2str(f1) 'Hz oraz ' num2str(f2) 'Hz'])
xlabel('Częstotliwość [Hz]'); ylabel('Amplituda [-]');
set(findobj(gcf,'type','axes'),'FontName','Halvetica','FontSize',12);