Zadanie 1
W tej części laboratorium należało stworzyć własne implementacje dyskretnej transformaty Fouriera oraz przetestować ich działanie dla sygnału typu “chirp”, o częstotliwości początkowej i końcowej dla czasu , przy częstotliwości próbkowania i całkowitym czasie .
Kod:
Fs = 1000; % Częstotliwość próbkowania
dt = 1/Fs; % Krok czasowy
T = 1; % Całkowity czas
t = 0:dt:T-dt; % Wektor czasu
f_start = 20;
f_stop = 220;
x = chirp(t,f_start,t(end),f_stop); % Sygnał chirp
N = length(x); % Ilość próbek
n = 0:(N-1); % Kolejne próbki
% Transformata Fouriera:
% Moja pierwsza implementacja:
for k=0:(N-1)
Core = cos(-2*pi*k*n/N)+1i*sin(-2*pi*k*n/N); % Rdzeń transformaty
X1(k+1) = sum(x.*Core);
end
% Moja druga implementacja:
for k=0:(N-1)
Core = exp(-2*pi*1i*k*n/N); % Rdzeń transformaty
X2(k+1) = sum(x.*Core);
end
% Funkcja Matlaba:
X3 = fft(x);
df = Fs/N; % Rozdzielczość spektralna
f_vec = (0:N-1)*df; % Wektor częstotliwości
X1_amp = abs(X1); % Amplituda
X1_phase = atan2(imag(X1),real(X1)); % Faza
X2_amp = abs(X2);
X2_phase = atan2(imag(X2),real(X2));
X3_amp = abs(X3);
X3_phase = atan2(imag(X3),real(X3));
% Normalizacja wektorów:
X1_amp_norm = X1_amp/N*2; X1_amp_norm(1)=X1_amp_norm(1)/2;
X1_phase_norm = X1_phase; X1_phase_norm(X1_amp_norm<0.01)=0;
X2_amp_norm = X2_amp/N*2; X2_amp_norm(1)=X2_amp_norm(1)/2;
X2_phase_norm = X2_phase; X2_phase_norm(X2_amp_norm<0.01)=0;
X3_amp_norm = X3_amp/N*2; X3_amp_norm(1)=X3_amp_norm(1)/2;
X3_phase_norm = X3_phase; X3_phase_norm(X3_amp_norm<0.01)=0;
figure(1)
subplot(3,2,1)
plot(f_vec, X1_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/2]);
title('Widmo amplitudowe dla pierszej implementacji FFT')
subplot(3,2,2)
plot(f_vec, X1_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
% display y-axis ticks in radians:
y_ticks = -1:0.25:1;
yticks(y_ticks*pi)
ticks_for_phase={};
for i=1:length(y_ticks)
ticks_for_phase{i}=[num2str(y_ticks(i)) ' \pi'];
end
yticklabels(ticks_for_phase)
xlim([0,Fs/2]);
title('Widmo fazowe dla pierszej implementacji FFT')
subplot(3,2,3)
plot(f_vec, X2_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/2]);
title('Widmo amplitudowe dla drugiej implementacji FFT')
subplot(3,2,4)
plot(f_vec, X2_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
% display y-axis ticks in radians:
y_ticks = -1:0.25:1;
yticks(y_ticks*pi)
ticks_for_phase={};
for i=1:length(y_ticks)
ticks_for_phase{i}=[num2str(y_ticks(i)) ' \pi'];
end
yticklabels(ticks_for_phase)
xlim([0,Fs/2]);
title('Widmo fazowe dla drugiej implementacji FFT')
subplot(3,2,5)
plot(f_vec, X3_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/2]);
title('Widmo amplitudowe przy użyciu funkcji Matlaba')
subplot(3,2,6)
plot(f_vec, X3_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
% display y-axis ticks in radians:
y_ticks = -1:0.25:1;
yticks(y_ticks*pi)
ticks_for_phase={};
for i=1:length(y_ticks)
ticks_for_phase{i}=[num2str(y_ticks(i)) ' \pi'];
end
yticklabels(ticks_for_phase)
xlim([0,Fs/2]);
title('Widmo fazowe przy użyciu funkcji Matlaba')
Wynik działania programu:
Jak widać, uzyskane widma są identyczne dla obydwu implementacji oraz przy użyciu funkcji Matlaba.
Zadanie 2
W tej części laboratorium należało wyznaczyć widmo amplitudowe z sumy dwóch sinusoid o częstotliwościach oraz , przy częstotliwości próbkowania , dla całkowitego czasu i .
Kod:
Fs = 1000; % Częstotliwość próbkowania
dt = 1/Fs; % Krok czasowy
T1 = 1; % Całkowity czas sygnału 1
t1 = 0:dt:T1-dt; % Wektor czasu sygnału 1
T2 = 10; % Całkowity czas sygnału 2
t2 = 0:dt:T2-dt; % Wektor czasu sygnału 2
f1 = 50;
f2 = 51;
x1 = sin(2*pi*f1*t1) + sin(2*pi*f2*t1); % Sygnał 1
x2 = sin(2*pi*f1*t2) + sin(2*pi*f2*t2); % Sygnał 2
N1 = length(x1); % Ilość próbek sygnału 1
N2 = length(x2); % Ilość próbek sygnału 2
% Transformata Fouriera:
X1 = fft(x1);
X2 = fft(x2);
df1 = Fs/N1; % Rozdzielczość spektralna 1
f_vec1 = (0:N1-1)*df1; % Wektor częstotliwości 1
X1_amp = abs(X1); % Amplituda 1
df2 = Fs/N2; % Rozdzielczość spektralna 2
f_vec2 = (0:N2-1)*df2; % Wektor częstotliwości 2
X2_amp = abs(X2); % Amplituda 2
% Normalizacja:
X1_amp_norm = abs(X1)/N1*2; X1_amp_norm(1)=X1_amp_norm(1)/2;
X2_amp_norm = abs(X2)/N2*2; X2_amp_norm(1)=X2_amp_norm(1)/2;
figure(1)
% Widmo amplitudowe sygnału 1:
subplot(2,1,1)
plot(f_vec1, X1_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/10]);
title('Widmo amplitudowe sygnału dla T = 1s')
% Widmo amplitudowe sygnału 2:
subplot(2,1,2)
plot(f_vec2, X2_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/10]);
title('Widmo amplitudowe sygnału dla T = 10s')
Wynik działania programu:
Jak widać, dla całkowitego czasu rozdzielczość spektralna przyjmuje zbyt dużą wartość, przez co nie da się odróżnić dwóch komponentów częstotliwościowych.
Zadanie 3
W tej części laboratorium należało stworzyć własną implementację krótko-okienkowej transformaty Fouriera (STFT) i przetestować jej działanie na przykładzie funkcji chirp z zadania 1.
Kod:
Fs = 1000; % Częstotliwość próbkowania
dt = 1/Fs; % Krok czasowy
T = 1; % Całkowity czas
t = 0:dt:T-dt; % Wektor czasu
f_start = 20;
f_stop = 220;
x = chirp(t,f_start,t(end),f_stop); % Sygnał chirp
N = length(x); % Ilość próbek
%WinTime = 0.1; % Długość okna w sekundach
WinTime = 0.25;
%Overlap = 50; % Nakładanie się, w procentach
Overlap = 0;
[T_VEC, F_VEC, STFT] = mySTFT(Fs,x,WinTime,Overlap);
figure(1)
imagesc(T_VEC, F_VEC, STFT);
axis('xy') % Odwrócenie osi XY
title('Własna implementacja STFT');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2]) % Widmo jednostronne
function [T_VEC, F_VEC, STFT] = mySTFT(Fs, x, WinTime, Overlap)
WIN = WinTime*Fs; % Długość okna w próbkach
Overlap = Overlap/100; % Parametr nakładania, w procentach
PRZES = WIN*Overlap; % Parametr przesuwania
if PRZES==0
PRZES=WIN;
end
AP = 1; % Aktualna próbka
index=1;
false=1;
true=0;
while false~=true
KP = AP+WIN-1; % Końcowa próbka
% Sprawdzenie, czy nie wychodzimy poza zakres:
if KP > length(x)
break;
end
seg = x(AP:KP); % Fragment sygnału
T_VEC(index) = (AP+(WIN/2)-1)/Fs; % Wektor czasu
window = hann(WIN)'; % Okno Hanninga
win_seg = seg.*window;
FFT = fft(win_seg); % FFT z fragmentu sygnału
STFT(index,:) = abs(FFT); % Macierz STFT
AP = AP+PRZES; % Nowa wartość AP
index = index+1;
end
STFT = STFT';
F_VEC = (0:(WIN-1))*(Fs/WIN); % Wektor częstotliwości
end
Wynik działania programu:
- Długość okna czasowego: 0,25 s, 0% nakładania:
- Długość okna czasowego: 0,1 s, 50% nakładania:
Zadanie 4
Kod:
Fs = 100; % Częstotliwość próbkowania
dt = 1/Fs; % Krok czasowy
T = 4; % Całkowity czas sygnału
t = 0:dt:T-dt; % Wektor czasu sygnału
WinTime1 = 0.1; % Długość okna 1 w sekundach
WinTime2 = 2; % Długość okna 2 w sekundach
Overlap = 50; % Nakładanie się w procentach
f1 = 20;
f2 = 22;
x = sin(2*pi*f1*t) + sin(2*pi*f2*t); % Sygnał
N = length(x); % Ilość próbek sygnału
df = Fs/N; % Rozdzielczość spektralna
f_vec = (0:N-1)*df; % Wektor częstotliwości
X = fft(x);
X_amp = abs(X); % Amplituda
% Normalizacja:
X_amp_norm = abs(X)/N*2; X_amp_norm(1)=X_amp_norm(1)/2;
[T_VEC1, F_VEC1, STFT1] = mySTFT(Fs,x,WinTime1,Overlap);
[T_VEC2, F_VEC2, STFT2] = mySTFT(Fs,x,WinTime2,Overlap);
figure(1)
plot(t,x)
xlabel('Czas [s]')
ylabel('Amplituda')
title('Sygnał w dziedzinie czasu')
figure(2)
plot(f_vec, X_amp_norm)
xlabel('Częstotliwość [Hz]')
ylabel('Amplituda')
xlim([0,Fs/2]);
title('Widmo amplitudowe')
figure(3)
imagesc(T_VEC1, F_VEC1, STFT1);
axis('xy')
title('Spektrogram dla okna o długości t = 0,1 s');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2])
figure(4)
imagesc(T_VEC2, F_VEC2, STFT2);
axis('xy')
title('Spektrogram dla okna o długości t = 2 s');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2])
Wynik działania programu:
Zadanie 5
Kod:
Fs=1000;
t=0:1/Fs:10;
x1=chirp(t,20,10,220); % Chirp
x2=zeros(size(t));
x2(5000)=1; % Impuls Diraca
x3=sin(10*2*pi*t)+sin(100*2*pi*t)+sin(250*2*pi*t); % Suma sinusoid
WinTime = 1; % Długość okna w sekundach
Overlap = 50; % Nakładanie się w procentach
[T_VEC1, F_VEC1, STFT1] = mySTFT(Fs,x1,WinTime,Overlap);
[T_VEC2, F_VEC2, STFT2] = mySTFT(Fs,x2,WinTime,Overlap);
[T_VEC3, F_VEC3, STFT3] = mySTFT(Fs,x3,WinTime,Overlap);
figure(1)
imagesc(T_VEC1, F_VEC1, STFT1);
axis('xy')
title('Spektrogram dla sygnału chirp');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2])
figure(2)
imagesc(T_VEC2, F_VEC2, STFT2);
axis('xy')
title('Spektrogram dla impulsu Diraca');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2])
figure(3)
imagesc(T_VEC3, F_VEC3, STFT3);
axis('xy')
title('Spektrogram dla sumy sinusoid');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
ylim([0,Fs/2])
figure(4)
cwt(x1,'amor',Fs)
title('Transformata falkowa dla sygnału chirp');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
figure(5)
cwt(x2,'amor',Fs)
title('Transformata falkowa dla impulsu Diraca');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
figure(6)
cwt(x3,'amor',Fs)
title('Transformata falkowa dla sumy sinusoid');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
Wynik działania programu:
