Analiza i identyfikacja sygnałów
Lab. 4:Wprowadzenie do przetwarzania sygnałów w dziedzinie częstotliwości.
Tomasz Osuch MECHATRONIKA 23
Ćw.1.
Wygenerować kilkusekundowy sygnał będący złożeniem trzech przebiegów
sinusoidalnych o częstotliwościach: 10, 80 i 120 Hz. Amplitudy sygnałów powinny
znajdować się w proporcjach 1:3:1 oraz fazy powinny być przesunięte o ok. 20o.
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 80;
f3 = 120;
fs = 480; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu
% Definicje przykładowych przebiegów sinusoidalnych
x1 = sin(2*pi*f1*t);
x4 = 3*sin(2*pi*f2*t+pi/9);
x6 = sin(2*pi*f3*t+2*(pi/9));
wynik=x1+x4+x6;
plot(t,wynik)
Ćw.2:
Utworzyć transformatę Fouriera sygnału wygenerowanego w p. 2 i pokazać ją graficznie w postaci wykresów amplitudowo-częstotliwościowych oraz fazowo-częstotliwościowych. Należy zwrócić uwagę na dobór odpowiedniej częstotliwości próbkowania oraz poprawne wyskalowanie osi częstotliwości. Na wykresie zaznaczyć pionową linią częstotliwość Nyqiusta. Skomentować uzyskane wyniki. (Uwaga: Do zestawiania wyników transformacji można wykorzystać funkcję „subplot”).
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 80;
f3 = 120;
fs = 40;
t=(0:1/fs:1);
n=fs/2;
% Definicje przykładowych przebiegów sinusoidalnych
x1 = sin(2*pi*f1*t);
x2 = 3*sin(2*pi*f2*t+pi/9);
x3 = sin(2*pi*f3*t+2*(pi/9));
x=x1+x2+x3;
y=fft(x);
m=abs(y);
c=(0:0.1:70);
b=(-400:1:0)
p=unwrap(angle(y));
f=(0:length(y)-1)'*fs/length(y);
subplot(1,2,1), plot(f,m,'black',n,c,'red')
title 'amplitudowa';
set(gca,'XTick',[10 20 30 40]);
subplot(1,2,2), plot(f,p*180/pi,'black',n,b,'red')
title 'fazowa';
set(gca,'XTick',[10 20 30 40]);
Otrzymany wykres amplitudowo częstotliwościowy jest symetryczny zdłuż lini Nyguista czyli polowy częstotliwości próbkowania ( w tym przypadku 20Hz). Amplituda przyjmuje ekstrema przy częstotliwości 10 i 30Hz. Wykres częstotlowściowo-fazowy nie jest symetryczny, a wartośc fazy maleje wraz z wzrostem częstotliwości.
Ćw.3:
Zmienić częstotliwość próbkowania sygnału z zadania 2 tak aby częstotliwość Nyquista była równa 100 Hz oraz wykonać zadanie 2 dodatkowo usuwając sinusoidę o częstotliwości 80 Hz. Skomentować wynik działania skryptu.
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 80;
f3 = 120;
fs = 200;
t=(0:1/fs:1);
n=fs/2;
% Definicje przykładowych przebiegów sinusoidalnych
x1 = sin(2*pi*f1*t);
x2 = 3*sin(2*pi*f2*t+pi/9);
x3 = sin(2*pi*f3*t+2*(pi/9));
x=x1+x2+x3;
y=fft(x);
m=abs(y);
c=(0:0.1:250);
b=(-100:1:500)
p=unwrap(angle(y));
f=(0:length(y)-1)'*fs/length(y);
subplot(1,2,1), plot(f,m,'black',n,c,'red')
title 'amplitudowa';
set(gca,'XTick',[50 100 150 200]);
subplot(1,2,2), plot(f,p*180/pi,'black',n,b,'red')
title 'fazowa';
set(gca,'XTick',[50 100 150 200]);
Pkt. 2.
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
%f2 = 80;
f3 = 120;
fs = 200;
t=(0:1/fs:1);
n=fs/2;
% Definicje przykładowych przebiegów sinusoidalnych
x1 = sin(2*pi*f1*t);
%x2 = 3*sin(2*pi*f2*t+pi/9);
x3 = sin(2*pi*f3*t+2*(pi/9));
x=x1+x3;
y=fft(x);
m=abs(y);
c=(0:0.1:100);
b=(-200:1:200)
p=unwrap(angle(y));
f=(0:length(y)-1)'*fs/length(y);
subplot(1,2,1), plot(f,m,'black',n,c,'red')
title 'amplitudowa';
set(gca,'XTick',[50 100 150 200]);
subplot(1,2,2), plot(f,p*180/pi,'black',n,b,'red')
title 'fazowa';
set(gca,'XTick',[50 100 150 200]);
Usuwając składową 80Hz otrzymalismy wykres symetryczny tak jak z tą składową. Wartości maksymalne i wysokie pojawiają się w przy tych samych częstotliwościach jednak ich wartości pod wpływem zmiany składowej zmieniły się. Natomiast wykresy fazowe różnią się znacznie. Dla pierwszej wartości sygnału prosta Nyquista przecina się z wykresem na wysokości 180 a dla drugiego na wysokości 0.
Ćw.4:
Dokonać odwrotnej transformacji Fouriera sygnału z zadania 2 oraz porównać wyniki ze sygnałem oryginalnym. Porównanie wykonać poprzez wykonanie obu wykresów w jednym oknie.
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 80;
f3 = 120;
fs = 200;
t=(0:1/fs:1);
NY=fs/2;
% Definicje przykładowych przebiegów sinusoidalnych
x1 = sin(2*pi*f1*t);
x2 = 3*sin(2*pi*f2*t+pi/9);
x3 = sin(2*pi*f3*t+2*(pi/9));
x=x1+x2+x3;
y=fft(x);
y1=ifft(x);
m=abs(y);
m1=abs(y1)
c=(0:0.1:250);
b=(-500:1:500)
p=unwrap(angle(y));
p1=unwrap(angle(y1));
f=(0:length(y)-1)'*fs/length(y);
f1=(0:length(y1)-1)'*fs/length(y1);
subplot(1,2,1), plot(f,m,'black',f1,m1,'red',NY,c,'green')
title 'amplitudowa';
set(gca,'XTick',[50 100 150 200]);
subplot(1,2,2), plot(f,p*180/pi,'black',f1,p1*180/pi,'red',NY,b,'green')
title 'fazowa';
set(gca,'XTick',[50 100 150 200]);
Jak widzimy po wykresie fazowo-częstotliwościowym te dwie funkcje są do siebie odwrotne.
Ćw.5:
Wykonać transformacje Fouriera przykładowych sygnałów zawartch w podpunkcie: „Przykład transformacji Fouriera dla wybranych sygnałów”. Dodatkowo wykonać następujące transformaty:
• dwóch impulsów występujących w bliskiej odległości
• przebiegu szumu losowego
• przebiegu kroku jednostkowego
• przebiegu trójkątnego
Ćw.6:
Wygenerować sygnał sinusoidalny o częstotliwości 11 Hz próbkowany z częstotliwością 200Hz. Przygotować dwie realizacje tego sygnału jedną o długości 1s oraz drugą o długości 0.95s. Wykonać przekształcenie Fouriera dla obu próbek przygotowanego sygnału. Moduły obu transformat narysować na jednym wykresie. Skomentować uzyskane wyniki.
f = 11;
fs = 200; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu
t1=0:(1/fs):0.95; % Wektor czasu
x = sin(2*pi*f*t);
x1 = sin(2*pi*f*t1);
y=fft(x);
y1=fft(x1);
m=abs(y);
m1=abs(y1)
p=unwrap(angle(y));
p1=unwrap(angle(y1));
f=(0:length(y)-1)'*fs/length(y);
f1=(0:length(y1)-1)'*fs/length(y1);
subplot(1,2,1), plot(f,m,'black',f1,m1,'red')
title 'amplitudowa';
set(gca,'XTick',[50 100 150 200]);
subplot(1,2,2), plot(f,p*180/pi,'black',f1,p1*180/pi,'red')
title 'fazowa';
set(gca,'XTick',[50 100 150 200]);
Dla mniejszego czasu sygnału wykres ma większą amplitude- nie zmienia się jednak częstotliwośc przy której te ekstrema amplitudy występuja. Przebieg fazowy dla większego czasu jest bardziej gładki jednak zawiera się w większym przedziale.
Ćw.7:
Powtórzyć eksperyment z punktu 6 nakładając wcześniej na sygnały okno Hanninga.
Skomentować uzyskane wyniki.
f = 11;
fs = 200; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu
t1=0:(1/fs):0.95; % Wektor czasu
x=sin(2*pi*f*t);
x1 = sin(2*pi*f*t1);
% Zastosowanie okna Hanninga
w = hanning(length(x));
w1 = hanning(length(x1));
y = fft(x'.*w);
y1 = fft(x1'.*w1);
f = (0:length(y)-1)'*fs/length(y);
f1 = (0:length(y1)-1)'*fs/length(y1);
ind=find(f<=fs/2);
ind=find(f1<=fs/2);
f2=f(ind); am2 = 20*log10(abs(y(ind)/norm(w)));
f3=f1(ind); am3 = 20*log10(abs(y1(ind)/norm(w1)));
plot(f2,am2,'black',f3,am3,'red')
title 'Hanninga'
Ćw.8:
Zmodyfikować sygnał testowy tworzony w ramach zadania 2 poprzez dodanie do niego sinusoid o częstotliwościach 12 i 81 Hz. Wykonać porównanie działania różnych okien czasowych nakładanych na testowy sygnał. Sprawdzić jaki wpływ na długość zastosowanego okna na wynikowe widmo.
% Zastosowanie okien czasowych
fs = 200; % Częstotliwość próbkowania sygnału
t=0:(1/fs):1; % Wektor czasu
% Utworzenie przebiegu sinusoidalnego
f1 = 10; % Częstotliwości przykładowych przebiegów sinusoidalnych
f2 = 80;
f3 = 120;
f4= 12;
f5= 81;
x1 = sin(2*pi*f1*t);
x2 = 3*sin(2*pi*f2*t+pi/9);
x3 = sin(2*pi*f3*t+2*(pi/9));
x4 = 3*sin(2*pi*f4*t+pi/9);
x5 = sin(2*pi*f5*t+2*(pi/9));
x=x1+x2+x3+x4+x5;
% Zastosowanie okna prostokątnego
y = fft(x);
f = (0:length(y)-1)'*fs/length(y);
ind=find(f<=fs/2);
f1=f(ind); am1 = 20*log10(abs(y(ind)));
figure(1)
plot(f1,am1)
title 'prostokątnego'
%pause;
% Zastosowanie okna Hanninga
w = hanning(length(x));
y = fft(x'.*w);
ind=find(f<=fs/2);
f2=f(ind); am2 = 20*log10(abs(y(ind)/norm(w)));
figure(2)
plot(f2,am2)
title 'Hanninga'
%pause;
% Zastosowanie okna Hamminga
w = hamming(length(x));
y = fft(x'.* w);
ind=find(f<=fs/2);
f3=f(ind); am3 = 20*log10(abs(y(ind)/norm(w)));
figure(3)
plot(f3,am3)
title 'Hamminga'
%pause;
% Zastosowanie okna Kaisera
w = kaiser(length(x),6);
y = fft(x'.* w);
ind=find(f<=fs/2);
f4=f(ind); am4 = 20*log10(abs(y(ind)/norm(w)));
figure(4)
plot(f4,am4)
title 'Kaisera'
%pause;
% Zastosowanie okna Chebyshewa
w = chebwin(length(x),40);
y = fft(x'.* w);
ind=find(f<=fs/2);
f5=f(ind); am5 = 20*log10(abs(y(ind)/norm(w)));
figure(5)
plot(f5,am5)
title 'Chebyshewa'
Wystąpiły poszarbania w mijescu maksymalnych amplitud.
Wystąpiło zniekształcenie pomiędzy wierzchołkami wykresu. W zamian za to wierzchołki wykresy są wygładzone.
Wyłagodzenie wirzchołków funkcji. Powstanie zniekształcenia przy podstawie ekstremów. Obcięcie łagodnych przejśc w wierzhołki funkcji.
Wyłagodzenie wirzchołków funkcji. Powstanie zniekształcenia przy podstawie ekstremów. Obcięcie łagodnych przejśc w wierzhołki funkcji.
Wyłagodzenie przebiegu funkjci.