Sprawozdanie - PSIWMUM, Lab4

Zadanie 1. Animacja wzoru Eulera

W tym zadaniu należało wykonać “animację” wzoru Eulera dla fazy w zakresie $< 0, 360° >$ , z krokiem $5°$ .

Kod:


faza=0:5:360;
y1=cos((faza./180).*pi); % Wykres cosinusa
y2=sin((faza./180).*pi); % Wykres sinusa

for tfaza=0:5:360 
    
    signal1=cos((tfaza/180)*pi); % Położenie punktu na wykresie cosinusa
    signal2=sin((tfaza/180)*pi); % Położenie punktu na wykresie sinusa

    figure(1)
    subplot(2,2,2);
    plot(faza,y1,'b-',tfaza,signal1,'ro');
    axis([0 360 -1 1]);
    ylabel('Amplituda [-]'), xlabel('Faza [°]'); 
    title('Cosinus');
    subplot(2,2,4);
    plot(faza,y2,'b-',tfaza,signal2,'ro');
    axis([0 360 -1 1]);
    ylabel('Amplituda [-]'), xlabel('Faza [°]'); 
    title('Sinus');
    subplot(2,2,[1,3]);
    compassplot((tfaza/180)*pi,1)
    title('Euler');
    pause(1)                                              
end

Wynik działania programu:

Zadanie 2. Własna implementacja DFT

W tym zadaniu należało napisać funkcję będącą własną implementacją Dyskretnej Transformaty Fouriera, przy użyciu jednej pętli for oraz funkcji sum().

Kod:


Fs=100;                             % Częstotliwość próbkowania
dt=1/Fs;                            % Krok czasowy
t=0:dt:10-1/Fs;                     % Wektor czasu
s=sin(2*pi*15*t) + sin(2*pi*40*t);  % Sygnał
N=length(s);                        % Ilość próbek
n=0:(N-1);                          % Kolejne próbki

% Moja implementacja Dyskretnej Transformaty Fouriera:
index=1;
for k=0:(N-1)
    S(index) = sum(s.*exp(-2*pi*1i*k*n/N));
    index=index+1;
end

f = (0:N-1)*Fs/N;                   % Wektor częstotliwości
S_amp = abs(S);                     % Amplituda
S_phase = atan2(imag(S),real(S));   % Faza

% Normalizuję wektory:
S_amp_norm = abs(S)/N*2; S_amp_norm(1)=S_amp_norm(1)/2;
S_phase_norm = S_phase; S_phase_norm(S_amp_norm<0.01)=0;

figure(1)
% Widmo amplitudowe:
subplot(2,1,1)
plot(f, S_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,50]);
% Widmo fazowe:
subplot(2,1,2)
plot(f, S_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,50]);

Wynik działania programu:

Na podstawie mojej implementacji Dyskretnej Transformaty Fouriera jestem w stanie utworzyć jednostronne widma amplitudowe i fazowe, z których można odczytać, że sygnał poddany DFT składa się z dwóch sygnałów o częstotliwościach odpowiednio 15 i 40 Hz, z przesunięciem fazowym $- \frac{π}{2}$ , co się zgadza, gdyż sygnał został zadany jako s=sin(2*pi*15*t) + sin(2*pi*40*t), a $cos (θ - \frac{π}{2}) = sin (θ)$ .

Zadanie 3. Widma częstotliwościowe sygnałów

W tym zadaniu należało przedstawić jednostronne widmo amplitudowe i widmo przesunięcia fazowego dla dwóch przypadków:

Przypadek 1: Zadany sygnał: $x (t) = 2 \cdot sin (2 π \cdot 26 \cdot t) + 4, 5 \cdot cos (2 π \cdot 80 \cdot t) + 0, 8 \cdot cos (2 π \cdot 8 \cdot t + π) + 6$

Kod:


Fs = 200;                           % Częstotliwość próbkowania
dt = 1/Fs;                          % Krok czasowy
t = 0:dt:3-dt;                      % Wektor czasu
% Sygnał:
s = 2*sin(2*pi*26*t) + 4.5*cos(2*pi*80*t)+0.8*cos(2*pi*8*t+pi)+6;
S = fft(s);                         % FFT dla sygnału
N = length(s);                      % Ilość próbek
df = Fs/N;                          % Krok częstotliwościowy
f_vec = (0:(N-1))*df;               % Wektor częstotliwości
S_amp = abs(S);                     % Amplituda
S_phase = atan2(imag(S),real(S));   % Faza

% Widma przed normalizacją:
figure(1)
% Widmo amplitudowe:
subplot(2,1,1)
plot(f_vec, S_amp)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,Fs/2]);
% Widmo fazowe:
subplot(2,1,2), plot(f_vec, S_phase)
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]);

% Normalizuję wektory:
S_amp_norm = abs(S)/N*2; S_amp_norm(1)=S_amp_norm(1)/2;
S_phase_norm = S_phase; S_phase_norm(S_amp_norm<0.01)=0;

% Widma po normalizacji:
figure(2)
% Widmo amplitudowe:
subplot(2,1,1) 
plot(f_vec, S_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/2]);
% Widmo fazowe:
subplot(2,1,2), plot(f_vec, S_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]);

Wynik działania programu:

Na podstawie uzyskanych widm można odczytać, że sygnał ma 4 składowe: $6 \cdot cos (2 π \cdot 0 \cdot t + 0 \cdot π)$ $1 \cdot cos (2 π \cdot 8 \cdot t + 1 \cdot π)$ $2 \cdot cos (2 π \cdot 26 \cdot t - \frac{1}{2} \cdot π)$ $4, 5 \cdot cos (2 π \cdot 80 \cdot t + 0 \cdot π)$ Co się zgadza, gdyż $cos (0) = 1$ oraz $cos (θ - \frac{π}{2}) = sin (θ)$ .

Przypadek 2:

Kod:


% Parametry sygnału:
A1 = 5;
A2 = 3;
A3 = 1.5;
f1=0;
f2=15;
f3=40;
phi1=0;
phi2=pi/4;
phi3=-pi/2;

Fs = 100;                           % Częstotliwość próbkowania
dt = 1/Fs;                          % Krok czasowy
t = 0:dt:3-dt;                      % Wektor czasu
% Sygnał:
s = A1*cos(2*pi*f1*t+phi1)+A2*cos(2*pi*f2*t+phi2)+A3*cos(2*pi*f3*t+phi3);
S = fft(s);                         % FFT dla sygnału
N = length(s);                      % Ilość próbek
df = Fs/N;                          % Krok częstotliwościowy
f_vec = (0:(N-1))*df;               % Wektor częstotliwości
S_amp = abs(S);                     % Amplituda
S_phase = atan2(imag(S),real(S));   % Faza

% Widma przed normalizacją:
figure(1)
% Widmo amplitudowe:
subplot(2,1,1)
plot(f_vec, S_amp)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,Fs/2]);
% Widmo fazowe:
subplot(2,1,2), plot(f_vec, S_phase)
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]);

% Normalizuję wektory:
S_amp_norm = abs(S)/N*2; S_amp_norm(1)=S_amp_norm(1)/2;
S_phase_norm = S_phase; S_phase_norm(S_amp_norm<0.01)=0;

% Widma po normalizacji:
figure(2)
% Widmo amplitudowe:
subplot(2,1,1) 
plot(f_vec, S_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda'),
xlim([0,Fs/2]);
% Widmo fazowe:
subplot(2,1,2), plot(f_vec, S_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]);

Wynik działania programu:

Zadanie 4. Aliasing

W tym zadaniu należało wykreślić widma amplitudowe dla tego samego sygnału, przy częstotliwościach próbkowania równym 1000 Hz, 100 Hz i 10 Hz.

Kod:


% Częstotliwości próbkowania sygnałów:
fs1=1000;
fs2=100;
fs3=10;

% Sygnał 1:
dt=1/fs1;
out = sim('Lab2_10_model2.slx');
s1=out.simout';
t1=out.tout';

% Zmniejszenie częstotliwości próbkowania sygnału 1:
s2=downsample(s1,10);
s3=downsample(s1,100);
% s2=decimate(s1,10);
% s3=decimate(s1,100);

% Dostosowanie wektorów czasu do utworzonych sygnałów:
t2 = downsample(t1, 10);
t3 = downsample(t1, 100);

% Ilości próbek sygnałów:
N1 = length(s1);
N2 = length(s2);
N3 = length(s3);

% Kroki częstotliwościowe:
df1 = fs1/N1;
df2 = fs2/N2;
df3 = fs3/N3;

% Wektory częstotliowści:
f_vec1 = (0:(N1-1))*df1;
f_vec2 = (0:(N2-1))*df2;
f_vec3 = (0:(N3-1))*df3;

% Obliczenie FFT dla sygnałów:
S1 = fft(s1);
S2 = fft(s2);
S3 = fft(s3);

% Amplitudy sygnałów:
S1_amp = abs(S1);
S2_amp = abs(S2);
S3_amp = abs(S3);

% Normalizowanie sygnałów:
S1_amp_norm = abs(S1)/N1*2; S1_amp_norm(1)=S1_amp_norm(1)/2;
S2_amp_norm = abs(S2)/N2*2; S2_amp_norm(1)=S2_amp_norm(1)/2;
S3_amp_norm = abs(S3)/N3*2; S3_amp_norm(1)=S3_amp_norm(1)/2;

figure(1)
% Przebieg czasowy:
subplot(2,1,1)
plot(t1,s1);
xlim([0 1]); 
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Sygnał spróbkowany, fs= 1000");
% Widmo amplitudowe:
subplot(2,1,2)
plot(f_vec1, S1_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs1/2]);

figure(2)
% Przebieg czasowy:
subplot(2,1,1)
plot(t2,s2);
xlim([0 1]); 
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Sygnał spróbkowany, fs= 100");
% Widmo amplitudowe:
subplot(2,1,2)
plot(f_vec2, S2_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs2/2]);

figure(3)
% Przebieg czasowy:
subplot(2,1,1)
plot(t3,s3);
xlim([0 1]); 
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Sygnał spróbkowany, fs= 10");
% Widmo amplitudowe:
subplot(2,1,2)
plot(f_vec3, S3_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs3/2]);

Wynik działania programu:

$f_{s} = 1000 Hz$
$f_{s} = 100 Hz$
$f_{s} = 10 Hz$

Jak widać, dla częstotliwości próbkowania równej 10 Hz nie jest spełnione kryterium Shannona-Nyquista. Dodatkowo, na podstawie widma amplitudowego można zauważyć, że oryginalny sygnał o składowych częstotliwościowych 6 Hz i 8 Hz jest błednie identyfikowany jako sygnał o składowych częstotliwościowych 2 Hz i 4 Hz, a więc zjawisko aliasingu wpływa także na analizę widmową sygnału.

Zadanie 5. Widma sygnałów charakterystycznych

W tym zadaniu należało wykreślić widma sygnałów:

Szum,
Impuls,
Sygnał piłokształtny.

Kod:


Tk = 25;    % Czas symulacji
fs = 1e4;   % Częstotliwość próbkowania
dt=1/fs;    % Krok czasowy

% Wygenerowanie sygnałów:
out = sim('gen_signal2.slx');
s1=out.signal1;
s2=out.signal2;
s3=out.signal3;
t=out.tout';

N = length(s1);         % Ilość próbek
df = fs/N;              % Krok częstotliwościowy
f_vec = (0:(N-1))*df;   % Wektor częstotliwści

% Obliczenie FFT dla sygnałów:
S1 = fft(s1);
S2 = fft(s2);
S3 = fft(s3);

% Amplitudy sygnałów:
S1_amp = abs(S1);
S2_amp = abs(S2);
S3_amp = abs(S3);

% Fazy sygnałów:
S1_phase = atan2(imag(S1),real(S1));
S2_phase = atan2(imag(S2),real(S2));
S3_phase = atan2(imag(S3),real(S3));

% Normalizowanie sygnałów:
S1_amp_norm = abs(S1)/N*2; S1_amp_norm(1)=S1_amp_norm(1)/2;
S1_phase_norm = S1_phase; S1_phase_norm(S1_amp_norm<0.01)=0;

S2_amp_norm = S2_amp;
S2_phase_norm = S2_phase; S2_phase_norm(S2_amp_norm<0.01)=0;

S3_amp_norm = abs(S3)/N*2; S3_amp_norm(1)=S3_amp_norm(1)/2;
S3_phase_norm = S3_phase; S3_phase_norm(S3_amp_norm<0.01)=0;

figure(1)
% Przebieg czasowy:
subplot(3,1,1)
plot(t,s1); 
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Szum");
% Widmo amplitudowe:
subplot(3,1,2)
plot(f_vec, S1_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs/2]);
% Widmo fazowe:
subplot(3,1,3)
plot(f_vec, S1_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
xlim([0,fs/2]);

figure(2)
% Przebieg czasowy:
subplot(3,1,1)
plot(t,s2);
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Impuls");
% Widmo amplitudowe:
subplot(3,1,2)
plot(f_vec, S2_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs/2]);
% Widmo fazowe:
subplot(3,1,3)
plot(f_vec, S2_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
xlim([0,fs/2]);

figure(3)
% Przebieg czasowy:
subplot(3,1,1)
plot(t,s3); 
xlabel("Czas [s]")
ylabel("Amplituda [-]")
title("Sygnał piłokształtny");
% Widmo amplitudowe:
subplot(3,1,2)
plot(f_vec, S3_amp_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Amplituda')
xlim([0,fs/2]);
% Widmo fazowe:
subplot(3,1,3)
plot(f_vec, S3_phase_norm)
xlabel('Częstotliwość [Hz]'), ylabel('Faza [rad]')
xlim([0,fs/2]);

Wynik działania programu:

Szum:
Impuls:
Sygnał piłokształtny:

Notatki

Explorer

Sprawozdanie - PSIWMUM, Lab4

Zadanie 1. Animacja wzoru Eulera

Zadanie 2. Własna implementacja DFT

Zadanie 3. Widma częstotliwościowe sygnałów

Zadanie 4. Aliasing

Zadanie 5. Widma sygnałów charakterystycznych

Graph View

Table of Contents