Untitled

clear all;
close all;
clc;

[syg, fp] = audioread('sygnal_gr2.wav');

N = length(syg); %ilosc probek
i =(0:N-1)';
%fp = N % częstotliwośc próbkowania
Tp = 1 / fp;

figure;
t = 0:Tp:(7.26 - Tp);
plot(t,syg);
xlim([0 7.26])
title('Przebieg czasowy');
ylabel('moduł');
xlabel('t[s]');
hold on;

X_fft = Tp * fft(syg,N);
X_gwm = abs(X_fft(1:N/2)).^2/N/Tp;
figure;
plot((0:N/2-1)/N*fp,10*log10(X_gwm));
title('GWM');
ylabel('dB');
xlabel('f [Hz]');
hold on;


syg_10 = syg(1:10000);
N_10 = length(syg_10); %ilosc probek
i =(0:N_10-1)';
%fp = N % częstotliwośc próbkowania
Tp = 1 / fp;

figure;
t = 0:Tp:(7.26 - Tp);
plot(t(1:10000),syg_10);
% xlim([0 7.26])
title('Przebieg czasowy');
ylabel('moduł');
xlabel('t[s]');
hold on;

X_fft = Tp * fft(syg_10,N_10);
X_gwm = abs(X_fft(1:N_10/2)).^2/N_10/Tp;
figure;
plot((0:N_10/2-1)/N_10*fp,10*log10(X_gwm));
title('GWM');
ylabel('dB');
xlabel('f [Hz]');
hold on;

syg_10 = syg(10001:20000);
N_10 = length(syg_10); %ilosc probek
i =(0:N_10-1)';
%fp = N % częstotliwośc próbkowania
Tp = 1 / fp;

figure;
t = 0:Tp:(7.26 - Tp);
plot(t(10001:20000),syg_10);
% xlim([0 7.26])
title('Przebieg czasowy');
ylabel('moduł');
xlabel('t[s]');
hold on;

X_fft = Tp * fft(syg_10,N_10);
X_gwm = abs(X_fft(1:N_10/2)).^2/N_10/Tp;
figure;
plot((0:N_10/2-1)/N_10*fp,10*log10(X_gwm));
title('GWM');
ylabel('dB');
xlabel('f [Hz]');
hold on;

figure;
spectrogram(syg,4000,3000,0:3999,8000);

% filtracja

r=1.2;
b=1.1;
f0 = 118;
a = -2 * cos(2*pi*f0/fp);
not_1= [1 a 1];
not_2=[r*r a*b 1];

[hn,wn]=freqz(not_1,not_2);
figure;
plot(wn*fp/2/pi,20*log10(abs(hn)))
title('filtr typu notch');
ylabel('moduł')
xlabel('f [Hz]');

figure
filteredSyg = filter(not_1,not_2,syg);
t = 0:Tp:(7.26 - Tp);
plot(t,filteredSyg);
title('Przebieg czasowy');
ylabel('moduł');
xlabel('t[s]');
xlim([0 7.26])

a = roots(not_2);
plot(a, 'o');
hold on;
plot(exp(-j*(0:1023)/1024*2*pi), 'r'), grid on;
title("Bieguny filtra na płaszyczźnie Nyquista");
figure;


figure
X_fft = Tp * fft(syg,N);
X_gwm = abs(X_fft(1:N/2)).^2/N/Tp;
figure;
plot((0:N/2-1)/N*fp,10*log10(X_gwm));
hold on;
X_fft = Tp * fft(filteredSyg,N);
X_gwm = abs(X_fft(1:N/2)).^2/N/Tp;
plot((0:N/2-1)/N*fp,10*log10(X_gwm));
title('GWM');
ylabel('dB');
xlabel('f [Hz]')
legend('sygnał wejściowy','przefiltrowany filtrem notch');

figure
impz(not_1,not_2);