Celem ćwiczenia jest przeprowadzenie analizy widmowej dyskretnych sygnałów okresowych przy zastosowaniu szybkiego przekształcenie Fouriera (FFT) oraz zbadaniu wpływu różnych czynników na dokładność przetwarzania. Do analizy i wizualizacji otrzymanych wyników wykorzystano funkcje programu MATLAB

Rezmer, J

Polish

Repozytorium Eny Politechnika Wrocławska

 1
 
Zakład Elektrotechniki Teoretycznej
 
 
ver.3 
 
 
 
 
 
 
ĆWICZENIE  III 
 
ANALIZA WIDMOWA SYGNAŁÓW DYSKRETNYCH 
 
 
 ? (00)   
 
Celem ćwiczenia jest przeprowadzenie analizy widmowej dyskretnych sygnałów okresowych 
przy zastosowaniu szybkiego przekształcenie Fouriera (FFT) oraz zbadaniu wpływu różnych 
czynników na dokładność przetwarzania. Do analizy i wizualizacji otrzymanych wyników 
wykorzystano funkcje programu MATLAB. 
 
Wykonując ćwiczenie zwróć uwagę na oznaczenia ikonami poszczególnych akapitów tej 
instrukcji. Mają one na celu usprawnienie wykonywania ćwiczeń. 
Oznaczenia   
 
 ?           przykłady  
 ?           opis funkcji MATLABA 
 ?           zadania do wykonania. 
 ? (15)   oznacza, iż aktualne zadanie powinno być wykonywane w czasie nie 
późniejszym niż 15 minuta ćwiczeń.  
 
 
Czas wykonania ćwiczenia wynosi 180 minut.  
 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 2
1.  WSTĘP 
 
FFT jest wykorzystywana do transformacji sygnałów zmiennych w czasie do 
dziedziny częstotliwości, czyli do wyznaczenia częstotliwościowej zawartości sygnału. 
Przykład 
 
W MATLABie algorytm dyskretnego przekształcenia Fouriera jest 
zrealizowany jako funkcja fft 
 
Funkcje Y = fft(x) oraz x = ifft(Y) implementują prostą i odwrotną 
dyskretną transformatę Fouriera zgodnie z zależnościami: 
( 1)( 1)
1
( ) ( )
N
n k
N
n
Y k x n W − −
=
= ∑       UWAGA na indeksowanie od 1 a nie od 0 
( 1)( 1)
1
1( ) ( )
N
n k
N
n
x n X k W
N
− − −
=
= ∑  
gdzie: N
j
N eW
π2−
=  
 ?            
Y = fft(x) daje w wyniku dyskretną transformatę Fouriera (DFT) 
wektora x 
Y = fft(x,N) daje w wyniku N-punktową DFT; jeśli długość wektora x 
jest mniejsza niż N, wektor x jest uzupełniany zerami do długości N, 
jeśli większa od N, to wektor x jest odpowiednio skracany. 
 ?   
Przykład 1 
FFT wektora x obliczamy następująco: 
» x = [4 3 7 -9 1 0 0 0]'; 
» y = fft(x) 
» y = 
» 6.0000 
» 11.4853 - 2.7574i 
» -2.0000 -12.0000i 
» -5.4853 +11.2426i 
» 18.0000 
» -5.4853 -11.2426i 
» -2.0000 +12.0000i 
» 11.4853 + 2.7574i 
 
Wektor x jest rzeczywisty natomiast y jest zespolony. W przykładzie pierwszy wyraz 
reprezentuje składową stałą, a piąty wyraz odpowiada częstotliwości Nyquista. Ostatnie trzy 
wyrazy odpowiadają częstotliwościom ujemnym i są odpowiednio sprzężone z 2, 3 i 4 
wyrazem.  
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 3
Najważniejsze informacje o sygnale poddanym transformacji zawarte są w widmie 
amplitudowym i fazowym DFT. W MATLABie otrzymujemy te charakterystyki przy pomocy 
funkcji abs i angle. 
 
Tworzymy przebieg składający się z jednej składowej sinusoidalnej. Najpierw tworzymy 
wektor czasu t i przy jego pomocy tworzymy sygnał x. 
 ?   
» t=0:1/100:10-1/100;%jaka jest częstotliwość próbkowania? 
» x=sin(2*pi*40*t); 
» y=fft(x); %obliczanie DFT ciągu x 
» m=abs(y); % widmo amplitudowe 
» p=unwrap(angle(y)); %widmo fazowe;  
Funkcja unwrap służy do usunięcia skoków fazy większych niż π. 
 
Tworzymy oś częstotliwości i rysujemy wykres amplitudy i fazy 
» f=(0:length(y)-1)'*100/length(y); %oś czestotliwości 
» subplot(2,1,1), plot(f,m); %wykres amplitudy 
» ylabel('modul amplitudy'), grid on 
» subplot(2,1,2), plot(f,p*180/pi); %faza 
» ylabel('faza [stopnie]'), grid on 
» xlabel('Czestotliwosc [Hz]') 
 
 
 
Wykres amplitudy jest symetryczny względem częstotliwości Nyquista 50 Hz. Całkowita 
informacja o sygnale zawarta jest w przedziale od 0 do 50 Hz (dlaczego ?). 
 ? (10)   
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 4
ĆWICZENIE 1 
 ?            
A. Przeprowadzić analizę amplitudy i fazy (DFT) dla 
sygnału o jednej składowej sinusoidalnej o różnych 
częstotliwościach [f], amplitudach [A] oraz 
częstotliwościach próbkowania [fp] przy akwizycji 
pełnej liczby okresów sygnału. 
Cw.1A-D 
gr f A fp 
1 50 6 150 
2 60 5 180 
3 70 4 210 
4 80 3 240 
5 90 2 270 
6 100 1 300 
 
B. Zbadać wpływ akwizycji niepełnej liczby okresów 
badanego sygnału na wyniki FFT.  
C. Zbadać wpływ położenia okna pomiarowego FFT.(wpływ przesunięcia czasowego). 
D. Zbadać wpływ na wyniki analizy DFT uzupełnienia zerami (zero-padding). 
 ?   
Przykład 2 
 
DFT stosuje się powszechnie do znajdowania składowych częstotliwościowych sygnału 
czasowego z szumem. Częstotliwość próbkowania wynosi 1000 Hz. Utworzono sygnał 
składający się z dwóch składowych o częstotliwościach 50 Hz i 120 Hz i szumu o średniej 
zerowej. 
» t=0:0.001:0.6; 
» x1=sin(2*pi*50*t)+sin(2*pi*120*t); 
» x=x1+2*randn(size(t)); 
» plot(x(1:50)) 
» title('sygnal z szumem');grid on 
» xlabel('czas [s]') 
 
 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 5
 
Patrząc na przebieg czasowy trudno powiedzieć, z jakich składowych widmowych składa się 
sygnał. Obliczamy 512 próbek DFT sygnału x : 
 ?   
» Y=fft(x,512); 
 
Widmo mocy, czyli miarę zawartości poszczególnych składowych sygnału w dziedzinie 
częstotliwości obliczamy z zależności: 
 
» Pyy=Y.*conj(Y)/512; 
 
Rysujemy wykres pierwszych 257 punktów na osi częstotliwości: 
 
» f=1000*(0:256)/512; 
» plot(f,Pyy(1:257));grid on 
» title('Widmo mocy y') 
» xlabel('f (Hz)') 
 
Wykres ten przedstawia zawartość częstotliwościową sygnału od częstotliwości 0 Hz do 
częstotliwości Nyquista 500 Hz.  
 
 
 
 
 ? (50)   
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 6
ĆWICZENIE 2 
 ?            
Cw.2 
gr f1 f2 K 
1 30 130 1,2 
2 60 100 1,3 
3 40 110 2,3 
4 20 90 2,1 
5 70 120 1,5 
6 80 140 1,7 
 
 Przeprowadzić analizę widmową sygnału będącego 
sumą 2 składowych sinusoidalnych o częstotliwościach 
[f1,f2], o amplitudach w zakresie od 1 do 3 przy poziomie 
szumu K. Szum gaussowski dodać do sygnału za pomocą 
funkcji randn, np: 
»x=sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+ 
K*randn(size(t)) 
 ?   
 
Przykład 3 
 
Jednym ze sposobów analizy widmowej sygnału jest obliczenie tzw. periodogramu 
zgodnie z zależnością: 
 
[ ] [ ]
[ ]
2
1
2
1
1
1
( )
N
jk
kj
N
k
w k x k e
Ne
w k
N
P
ω
ω
−
=
=
=
∑
∑
; 
gdzie x[k] – próbki sygnału, w[k] – dyskretne wartości okna modyfikującego próbki  
Periodogram przedstawia ciągły rozkład gęstości mocy sygnału w dziedzinie częstotliwości ( 
power spectral density - psd) i może być wyznaczany z transformaty Fouriera. 
 
W Matlabie do obliczeń periodogramu można wykorzystać podstawową funkcję FFT, lub 
bardziej złożone PSD lub PERIODOGRAM 
 
[Pxx,f] = periodogram(x,window,nfft,fs) 
wyznacza PSD sygnału x. Domyślnie wektor x jest modyfikowany oknem 
prostokątnym ( window = BOXCAR) o tej samej długości co wektor x, nfft – 
szerokość okna pomiarowego fft, fs – częstotliwość próbkowania.  
 ?   
fs=1000;           % częstotliwość próbkowania 1000 Hz 
t=(0:fs)/fs;       % próbki o czasie trwania 1 s 
A=[1 2];           % amplitudy sinusoid 
f=[140;150];       % częstotliwości sinusoid 
x=A*sin(2*pi*f*t)+0.1*randn(size(t)); %sygnał 
[Pxx,f]=periodogram(x,[],1024,fs); %periodogram 
semilogy(f,Pxx); grid; 
 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 7
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
-80
-70
-60
-50
-40
-30
-20
-10
0
10
Czestotliw osc (kHz)
M
oc
/c
ze
st
. (
dB
/H
z)
Periodogram
 
Przeciek widma 
 
 Jedną z wad periodogramu jest nakładanie się widma (przeciek). Występuje on 
najwyraźniej w przypadku badania krótkich sekwencji danych. Jeśli weźmiemy tylko 100 
próbek sygnału x to otrzymamy:  
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
-80
-70
-60
-50
-40
-30
-20
-10
0
czest. (kHz)
m
oc
/c
ze
st
 (d
B/
H
z)
Periodogram
 
Rozdzielczość 
Rozdzielczość jest zdolnością do rozróżnienia dwóch blisko siebie położonych 
składowych częstotliwościowych. Warunkiem przybliżonym rozróżnienia dwóch składowych 
jest  
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 8
 
1 2
sff f f
N
Δ = − > . 
 
W przypadku nie spełnienia tego warunku (61 próbek), otrzymano: 
 
periodogram(x(1:61),[],1024,fs); %periodogram 
 
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
-80
-70
-60
-50
-40
-30
-20
-10
0
czest. (kHz)
m
oc
/c
ze
st
. (
dB
/H
z)
Periodogram
 
 
 ? (90)   
ĆWICZENIE 3 
 ?            
Sprawdzić wpływ częstotliwości próbkowania, długości sekwencji próbek sygnału na 
zjawisko przecieku widma i zdolności rozdzielczej periodogramu. Zbadać wpływ poziomu 
szumu na właściwości periodogramu. 
 ?   
Przykład 4 
 
Modyfikacja periodogramu polega na zastosowaniu czasowej funkcji okna do próbek sygnału 
w celu zmniejszenia niektórych niekorzystnych właściwości. 
Wyznaczono periodogram sygnału o długości 100 próbek 
 
[Pxx,f]=periodogram(x,blackman(length(x)),1024,fs); %periodogram 
  
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 9
 
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
-80
-70
-60
-50
-40
-30
-20
-10
0
czest. (kHz)
m
oc
/c
ze
st
. (
dB
/H
z)
Periodogram
 
Dostępne funkcje okna: 
bartlett, blackman , chebwin, hamming, hann, kaiser, triang 
 ? (120)   
ĆWICZENIE 4 
 ?            
Cw.4 
gr Okno 1 Okno2 
1 Bartlett Chebwin 
2 Hamming Kaiser 
3 Hann0 Tukeywin
4 Bartlett  Kaiser 
5 Hamming Tukeywin
6 Hann Chebwin 
 
Zbadać wpływ dwóch wybranych funkcji okna na 
właściwości periodogramu (przeciek, rozdzielczość). 
Badania przeprowadzić na sygnale zawierającym składowe 
o bliskich częstotliwościach oraz szum.  
 
 
 
 
 ?   
Przykład 5 
Analiza harmoniczna sygnału, wyznaczenie mocy harmonicznych i rekonstrukcja 
sygnału.  
 Rozłożono sygnał trójkątny o okresie T=5s i amplitudzie 1 na składowe harmoniczne. 
512 próbek sygnału zbadano w celu określenia jaka część mocy sygnału zawarta jest w 
kolejnych składowych harmonicznych. 
 
 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 10
» T=5; 
» N=512; 
» t=linspace(0,T,N+1);t=t(1:N); oś czasu 
» x=sawtooth(0.4*pi*t,0.5);sygnał trójkątny 
» X=fft(x,512);  
» PSD=X.*conj(X)/N; 
» [sum(PSD) norm(x)^2]  Dlaczego te wyniki są sobie równe? 
» ans = 
» 170.6719 
»[moc,ind]=sort(PSD); dzięki operacji sortowania można wyznaczyć 
składowe widma o największej mocy. 
» m=6; ind(N:-1:N-m+1)  identyfikacja 6 składowych o największej mocy 
» ans = 
»       512   2   510   4   508   6 
 
Rekonstrukcja sygnału z jego harmonicznych 
 
» X6=zeros(1,N); zarezerwowanie przestrzeni dla rekonstruowanego widma 
» h=[512 2 510 4 508 6]; 
» X6(h)=X(h); kopiowanie 6 składowych do rekonstruowanego widma sygnału 
» xr=ifft(X6); odwrotna FFT  
» plot(t,[x;xr]);grid on; wykres 
 
 
 
 
» proc=100*(norm(X6)/norm(X))^2;  
» proc = 
»    99.9276 
 
Wniosek: 99,9276 % mocy sygnału zawarte jest w jego sześciu harmonicznych. Są to 
harmoniczne o indeksach 2, 4 i 6 odpowiadające częstotliwościom 0,2 Hz 0,4 i 0,6 Hz oraz 
składowe sprzężone. 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 
 11
 
 ? (160)   
ĆWICZENIE 5 
 ?            
Przeprowadzić analizę i rekonstrukcję dowolnego sygnału prostokątnego (funkcja 
square) i piłokształtnego (funkcja sawtooth). Określić zależności pomiędzy właściwościami 
sygnału a liczbą harmonicznych koniecznych do ich możliwie dokładnej rekonstrukcji 
(syntezy).  
 
 ? (180)   
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Literatura 
 
Oppenheim, R.W. Schaffer: ‚„Cyfrowe przetwarzanie sygnałów”. WKŁ, Warszawa 1979. 
W.Borodziewicz, K.Jaszczak: „Cyfrowe przetwarzanie sygnałów”. WNT, Warszawa 1987. 
A. Biran, M. Breiner: Matlab 5 for Engineers, Addison-Wesley, 1999 
„Matlab User’s Guide” The MathWorks Inc.1996  
 
? Przykład                            ? Help Matlaba                                   ? Zadanie                                     ? (Czas) 


Analiza widmowa sygnałów dyskretnych

http://zet10.ipee.pwr.wroc.pl/getfile.py?docid=144&name=3_fft_v3&format=pdf&version=1

Analiza widmowa sygnałów dyskretnych

Abstract

Similar works

Full text

Available Versions

Repozytorium Eny Politechnika Wrocławska