%% Cyfrowe przetwarzanie sygnałów drgań
% 1) Zdefiniowanie parametrów rejestracji i wczytanie danych

clear all

fp=48000;                       % częstotliwość próbkowania
tk=5;                           % czas pomiaru
load dane_lab_3
z=yt(:,2);

% 2) Zdefiniowanie parametrów filtracji dolnoprzepustowej
% Tę część należy modyfikować zmieniając watości N, Fp, Rp, Rst

N   = 100;                      % rząd filtra
Fp  = 20000;                    % częstotliwość "krawędzi" pasma przepustowego
Fs  = fp;                       % częstotliwość próbkowania
Rp  = 0.01;                     % tętnienie międzyszczytowe sygnału w paśmie przepustowym w [dB]
Rst = 1e-4;                     % tłumienie sygnału w paśmie zaporowym, np. 0.0001 odpowiada tłumieniu na poziomie 80 dB

% 3) Charakterystyka filtra dolnoprzepustowego

eqnum = firceqrip(N,Fp/(Fs/2),[Rp Rst],'passedge');     % eqnum = wektor współczynników filtra
fvtool(eqnum,'Fs',Fs,'Color','White')                   % wizualizacja charakterystyki filtra

% 4) Filtracja dolnoprzepustowa

lowpassFilt = dsp.LowpassFilter('DesignForMinimumOrder',false, ...
    'FilterOrder',N,'PassbandFrequency',Fp,'SampleRate',Fs,...
    'PassbandRipple',Rp, 'StopbandAttenuation',80);
z_filt_d = lowpassFilt(z);

% 5) Wygenerowanie wykresów drgań

figure
plot(yt(:,1),z,'b');
xlabel('czas[s]'); ylabel('y'); hold on
plot(yt(:,1),z_filt_d,'r');
legend('surowe dane','fitracja dolnoprzepustowa')

% 6) Widmo drgań

figure
Z=abs(fft(z))/((tk*fp)/2);                      % szybka transformata Fouriera sygnału drgań
Z_filt_d=abs(fft(z_filt_d))/((tk*fp)/2);        % szybka transformata Fouriera sygnału drgań po filtracji dolnoprzepustowej
f=0:1/tk:fp/2-1/tk;                             % wygenerowanie osi częstotliwości
plot(f,Z(1:(tk*fp)/2),'b');hold on;             % wykres widma sygnału drgań
plot(f,Z_filt_d(1:(tk*fp)/2),'r');hold on;      % wykres widma sygnału drgań po filtracji dolnoprzepustowej
xlabel('f [Hz]')                                % etykieta osi poziomej
ylabel('Z [m/s^2]')                             % etykieta osi pionowej
grid;                                           % siatka
xlim([0, fp/2])                                 % zakres wyswietlania danych na osi x
legend('surowe dane','fitracja dolnoprzepustowa')

% 7) Zdefiniowanie parametrów filtracji górnoprzepustowej
% Tę część należy modyfikować zmieniając watości Fpass, Fstop, Rp, Astop

filtertype = 'IIR';             % typ filtra
	Fpass = 2;                  % początek pasma przepustowego
	Fstop = 1;                  % granica pasma zaporowego
	Rp = 0.1;                   % tętnienie międzyszczytowe sygnału w paśmie przepustowym w [dB] 0.1 odpowiada 0.01 w [dB]
	Astop = 80;                 % tłumienie sygnału w paśmie zaporowym w [dB]

% 8) Charakterystyka filtra górnoprzepustowego

highpassFilt = dsp.HighpassFilter('SampleRate',Fs,...
                             'FilterType',filtertype,...
                             'PassbandFrequency',Fpass,...
                             'StopbandFrequency',Fstop,...
                             'PassbandRipple',Rp,...
                             'StopbandAttenuation',Astop);
fvtool(highpassFilt)                                    % wizualizacja charakterystyki filtra

% 9) Filtracja górnoprzepustowa

z_filt_dg = highpassFilt(z_filt_d);

% 10) Wygenerowanie wykresów drgań

figure
plot(yt(:,1),z,'b');
xlabel('czas[s]'); ylabel('y'); hold on
plot(yt(:,1),z_filt_d,'r');
plot(yt(:,1),z_filt_dg,'g');
legend('surowe dane','fitracja dolnoprzepustowa','fitracja dolno- i górno- przepustowa')

% 11) Widmo drgań

figure
Z=abs(fft(z))/((tk*fp)/2);                      % szybka transformata Fouriera sygnału drgań
Z_filt_d=abs(fft(z_filt_d))/((tk*fp)/2);        % szybka transformata Fouriera sygnału drgań po filtracji dolnoprzepustowej
Z_filt_dg=abs(fft(z_filt_dg))/((tk*fp)/2);      % szybka transformata Fouriera sygnału drgań po filtracji górno- i dolno- przepustowej
f=0:1/tk:fp/2-1/tk;                             % wygenerowanie osi częstotliwości
plot(f,Z(1:(tk*fp)/2),'b');hold on;             % wykres widma sygnału drgań
plot(f,Z_filt_d(1:(tk*fp)/2),'r');hold on;      % wykres widma sygnału drgań po filtracji dolnoprzepustowej 
plot(f,Z_filt_dg(1:(tk*fp)/2),'g');hold on;      % wykres widma sygnału drgań po filtracji górno- i dolno- przepustowa
xlabel('f [Hz]')                                % etykieta osi poziomej
ylabel('Z [m/s^2]')                             % etykieta osi pionowej
grid;                                           % siatka
xlim([0, fp/2])                                 % zakres wyswietlania danych na osi x
legend('surowe dane','fitracja dolnoprzepustowa','fitracja dolno- i górno- przepustowa')