zad1.m:
x=[1 0 0 0];
h=[0 2 1
0];
y=splot(h,x);
subplot(3,1,1)
stem(x)
subplot(3,1,2)
stem(h)
subplot(3,1,3)
stem(y)
splot.m:
function y=splot(h,x)
N=length(x);
M=length(h);
y=zeros(M+N-1,1);
for n=1:M+N-1
disp(sprintf('y(%d)=',n))
ind_x =max(1,n-M+1):min(N,n)
ind_h =min(n,M):-1:max(1,n-N+1)
y(n)=sum(h(ind_h).*x(ind_x));
end
x=[0 1 0 0 0]; h=[0 2 1 0 0];
x=[0 0 1 0 0]; h=[0 2 1 0 0];
x=[0 2 -1 0 0]; h=[0 -1 2 1 0];
przemienność
x(n)*y(n)=y(n)*x(n)
łączność
[x(n)*y(n)]*z(n)=x(n)*[y(n)*z(n)]
rozdzielność względem dodawania
x(n)*[y(n)+w(n)]=x(n)*y(n)+x(n)*w(n)
Przygotować funkcje zwracającą sgnał o zadanej długości w czasie i zadanej częstości próbkowania.
Potrzebne wersje funkcji to:
sin o zadanej częstości (w Hz) i fazie w stopniach
gabor ( gauss modulowany cosinusem) o zadanej częstości i standardowym odchyleniu w czasie
biały szum o zadanym odchyleniu standardowym
Wygenerować sygnał s próbkowany 128 Hz złożony z 2 sinusów o częstościach 10 i 30 Hz.
Obliczyć jego transformatę fft(s). Obejrzeć wynik.
Jaka jest interpretacja poszczególnych pików?
Dla sygnału s wyskalować czas.
Obliczyć i wyświetlić widmo sygnału s przy pomocy fft, wyskalować częstość.
Z transformaty odtworzyć oryginalny sygnał.
Z transformaty usunąć jeden z pików i odtworzyć sygnał.
Obejrzeć transformatę okna prostokątnego, Hanninga, kaisera
Prosze zrobic okno prostokatne (dlugosc 512 probek, jedynki od 200 do 300)
Okno transformujemy do dziedziny czestosci, w transformacie wydzielamy wartosc bezwzgledna (abs) i faze (funkcja angle)
wyswietlamy wartosci bezwzgledne i faze (zapoznac sie z funkcja unwrap)
przesuwamy faze o stala wartosc na dwa psosoby raz po unwapie a raz bez unwrapa
skladamy zaburzone transformaty spowrotem do liczb zespolonych (funkcja complex)
transformujemy do dziedziny czasu i ogladamy wynik. Co ciekawego widac? Prosze sprawdzic co robia jakies inne nieliniowe zaburzenia fazy
Maksymalna częstość w sygnale, a częstość próbkowania.
Ustalić częstość próbkowania na 32Hz
Wytworzyć sinusy o częstościach f od 2 do 32 Hz co 4Hz. Narysujcie sygnały i ich widma . Czy zaobserwowaliście jakościowe zmiany?
Jaka jest interpretacja poszczególnych pików?
Jakie stąd płyną wnioski co do próbkowania sygnałów ciągłych i przepróbkowywania dyskretnych z większej częstści prókowania na mniejszą?
Bardzo często musimy oszacować widmo mocy sygnału zawierającego znaczny udział szumu.
Poniższy skrypt ilustruje niepewność szacowania pików w widmie otrzymanym z transformaty fouriera dla sygnału zawierającego szum.
t=0:1/256:20-1/256;
s=sin(2*pi*t*30);
sz=randn(1,20*256);
s_sz=s+sz;
for i=1:20
x=s_sz(1+(i-1)*256:i*256);
X=fft(x);
sp=X.*conj(X)/256;
Psp(:,i)=sp(1:128);
end
plot(Psp)
Proszę obejrzeć otrzymane widma.
Jaka jest niepewność wniku?
Czy podobny problem występuje dla sygnału bez szumu?
Skonstruuj funkcję rysującą średnie widmo wraz z 95% przedziałem ufności.
Prosze wczytac sygnal test_signal
korzystajac z przygotowanej w domu funkcji narysowac widmo mocy -> okreslic czestosci wystepujace w sygnale i jego energie
Zrobić odcinek czasu t= [0 2) s próbkowany 128 Hz
Wygenerowac (s1 = ) sinusa f= 20 Hz okienkowanego gausem std =0.1 i pozycji t=1.5 s
Wygenerowac (s2 = ) sinusa f=20 Hz okienkowanego gausem std =0.1 i pozycji t = 0.5 s
Zapoznac sie z funkcja xcorr. Obejrzec funkcje autokorelacji s1 i s2
Zmienic t od 0.5:0.1:1 dla sygnalu s2 i dla kazdego z tych t obejrzec funkcje korelacji wzajemnej s1 i s2
Zmieniac f od 20Hz do 50Hz co 10Hz w obu sygnalach i obserwowac zachowanie funkcji korelacji wzajemnej
1) funkcja autokorelacji sygnalu okresowego z szumem zawiera informacje o jego okresowości
dt=1/256
t=0:dt:20-dt;
s=sin(2*pi*t*30);
sz=1.5*randn(1,20*256);
s_sz=s+sz;
for i=1:20
x=s_sz(1+(i-1)*256:i*256);
subplot(211)
plot(x)
[ak, lags] =xcorr(x);
subplot(212)
plot(lags*dt,ak)
title(i)
pause(1)
end
2) Sprawdzenie twierdzenia Wienera-Chińczyna: transformata Fouriera funkcji autokorelacji jest równa kwadratowi modułu transformaty Fouriera
fs=128; %czestosc probkowania [Hz]
dt = 1/fs; %przedzial probkowania [s]
fN=fs/2; %czestosc Nyquista
tmin = 0;
tmax = 2-dt;
t=tmin:dt:tmax; %[s]
sigma = 0.1;
t1 = 1.5;
okno1 = normpdf(t,t1,sigma);
f = 20;
s=sin(2*pi*t*f);
s1 = okno1.*s;
[ak, lags] = xcorr(s1,s1);
figure;
subplot(211)
plot(t,s1)
title('Sygnal s1 ');
subplot(212)
plot(lags*dt,ak);
title('Funkcja autokorelacji sygnalu s1 ');
xlabel('t[s]');
pause
% 9. Sprawdzenie twierdzenia Wienera - Chinczyna :
trans1 = fft (s1);
powerF= trans1 .* conj(trans1);
moduls1 = powerF(1:end/2+1);
transak = fft(ak);
% Transformate poddaje normalizacji i obliczam jej modul :
transak = transak(1:end/2+1);
transak = transak;
transak = (transak.*conj(transak)).^0.5
figure(9);
subplot(2,1,1);
stem((0:length(transak)-1)/length(transak)*fN,transak) ;
title('Modul transformaty funkcji autokorelacji');
xlabel('[Hz]');
subplot(2,1,2);
stem((0:length(moduls1)-1)/length(moduls1)*fN,moduls1) ;
title('Kwadrat transformaty sygnalu');
xlabel('[Hz]');
%end;
przykład sygnału wytworzonego z modelu AR:
a(1)=0.7;
a(2)=-0.3;
N=256;
x=zeros(1,N);
for i=3:N
x(i)=a(1)*x(i-1)+a(2)*x(i-2)+randn();
end
% estymacja parametrów modelu
subplot(211)
plot(x)
subplot(212)
ak=xcorr(x,'coeff');
plot(ak)
R=ak(N:end);
r0=R(1);
r1=R(2);
r2=R(3);
a2=(r1^2-r0*r2)/(r1^2-r0);
a1=r2/r1-r0/r1*a2;
disp([ a(1) a(2)])
disp([a1 a2])
widmo modelu:
Model:
lub zbierając wyrazy z
sygnałem po lewej stronie:
Po transformacji do dziedziny częstości:
zatem:
widmo mocy to
czyli że można je policzyć wprost ze współczynników modelu.