Krystian Babioch AiR sem. 5

P odstawy Cyfrowego Przetwarzania Sygnałów

Uzyskane wykresy:

M -file:

SygnalEKG = importdata('ekg0.mat'); %Import danych z pliku

[WymiarX, WymiarY] = size(SygnalEKG); %Odczyt wymiarów macierzy

CzestProb = 250; %Deklaracja czestotliwosci probkowania

OkresProb=1/CzestProb; %Obliczenie okresu probkowania

Czas=(0:OkresProb:WymiarX/CzestProb-OkresProb)'; %Obliczenie wartosci wektora czasu

%Zadanie 1.

figure; %Deklaracja nowego obszaru kreslenia

plot(Czas, SygnalEKG); %Rysowanie przebiegu czasowego podanego sygnału

axis([-Inf ,5 ,-Inf ,Inf]); %Ustawienie limitów osi

grid on; %Włączenie siatki na wykresie

title('Przebieg czasowy oryginalnego sygnału EKG') %Nadanie nazwy wykresu

legend('Sygnal EKG'); %Ustawienie etykiet legendy

SygnalEKGFourier=fft(SygnalEKG, CzestProb); %Dokonanie transformaty Fouriera na sygnale

SygnalGWM=abs(SygnalEKGFourier).^2./CzestProb; %Obliczenie wartości GWM

CzestNyquista = CzestProb/2; %Obliczenie częstotliwości Nyquista

figure; %Deklaracja nowego obszaru kreslenia

semilogy((0:CzestNyquista-1), SygnalGWM(1:CzestNyquista)); %Rysowanie gęstości widmowej mocy sygnału EKG

grid on; %Włączenie siatki na wykresie

title('Gęstość widmowa mocy oryginalnego sygnału EKG') %Nadanie nazwy wykresu

legend('Sygnal EKG'); %Ustawienie etykiet legendy

%Zadanie 2.

CzestZakl = 50; %Częstotliwość zakłócenia

CzestProbZakl = 250; %Częstotliwości z jaką jest próbkowane zakłócenie

OkresZakl = 1/CzestProbZakl; %Okres próbkowania zakłócenia

w = 2*pi*CzestZakl; %Obliczenie omegi zakłócenia

FiltrZerujacyLicznik = [1 -2*cos(w*OkresZakl) 1]; %Współczynniki licznika filtra zerującego

FiltrZerujacyMianownik = 1; %Współczynniki mianownika filtra zerującego

fvtool(FiltrZerujacyLicznik, FiltrZerujacyMianownik); %Rysowanie charakterystyk filtra zerującego

SygnalEKGPoFiltrZer = filter(FiltrZerujacyLicznik, FiltrZerujacyMianownik, SygnalEKG); %Filtracja sygnału filtrem zerującym

figure; %Deklaracja nowego obszaru kreslenia

plot(Czas, SygnalEKG, 'g'); %Rysowanie przebiegu czasowego oryginalnego sygnału

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, SygnalEKGPoFiltrZer, 'r'); %Rysowanie przebiegu czasowego przefiltrowanego sygnału

axis([-Inf ,5 ,-Inf ,Inf]); %Ustawienie limitów osi

grid on; %Włączenie siatki na wykresie

title('Prównanie przebiegó czasowych po przefiltrowaniu filtrem zerującym oraz przed')%Nadanie nazwy wykresu

legend('Sygnal EKG', 'Filtr Zerujacy'); %Ustawienie etykiet legendy

SygnalEKGPoFiltrZerFourier=fft(SygnalEKGPoFiltrZer,CzestProb); %Transformata Fouriera na przefiltrowanym sygnale

PoZerGWM=abs(SygnalEKGPoFiltrZerFourier).^2./CzestProb; %Obliczenie wartości GWM dla przefiltrowanego sygnału

CzestNyqZakl = CzestProbZakl/2; %Obliczenie częstotliwości Nyquista zakłócenia

figure; %Deklaracja nowego obszaru kreslenia

semilogy((0:CzestNyquista-1), SygnalGWM(1:CzestNyquista), 'Color', 'g'); %Rysowanie gęstości widmowej mocy sygnału EKG

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

semilogy((0:CzestNyqZakl-1), PoZerGWM(1:CzestNyqZakl), 'Color', 'r'); %Rysowanie gęstości widmowej mocy sygnału EKG po filtracji filtrem zerującym

grid on; %Włączenie siatki na wykresie

title('Prównanie gęstości widmowej mocy po przefiltrowaniu filtrem zerującym oraz przed')%Nadanie nazwy wykresu

legend('Sygnal EKG', 'Filtr Zerujacy'); %Ustawienie etykiet legendy

%Zadanie 4.

a1=0.7; %Dobierany współczynnik 1

a2=0.5; %Dobierany współczynnik 2

a3=0.9; %Dobierany współczynnik 3

FiltrWyglLicznik1 = 1-a1; %Współczynniki licznika filtra zerującego 1

FiltrWyglMianownik1 = [1 -a1]; %Współczynniki mianownika filtra zerującego 1

FiltrWyglLicznik2 = 1-a2; %Współczynniki licznika filtra zerującego 2

FiltrWyglMianownik2 = [1 -a2]; %Współczynniki mianownika filtra zerującego 2

FiltrWyglLicznik3 = 1-a3; %Współczynniki licznika filtra zerującego 3

FiltrWyglMianownik3 = [1 -a3]; %Współczynniki mianownika filtra zerującego 3

PoWygl1 = filter(FiltrWyglLicznik1, FiltrWyglMianownik1, SygnalEKGPoFiltrZer);%Sygnał po wygładzeniu podanym filtrem 1

PoWygl2 = filter(FiltrWyglLicznik2, FiltrWyglMianownik2, SygnalEKGPoFiltrZer);%Sygnał po wygładzeniu podanym filtrem 2

PoWygl3 = filter(FiltrWyglLicznik3, FiltrWyglMianownik3, SygnalEKGPoFiltrZer);%Sygnał po wygładzeniu podanym filtrem 3

figure; %Deklaracja nowego obszaru kreslenia

plot(Czas, SygnalEKGPoFiltrZer, 'y'); %Rysowanie przebiegu czasowego przefiltrowanego sygnału

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoWygl1, 'g'); %Rysowanie przebiegu czasowego wygładzonego sygnału 1

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoWygl2, 'r'); %Rysowanie przebiegu czasowego wygładzonego sygnału 2

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoWygl3, 'b'); %Rysowanie przebiegu czasowego wygładzonego sygnału 3

axis([-Inf ,5 ,-Inf ,Inf]); %Ustawienie limitów osi

grid on; %Włączenie siatki na wykresie

title('Porównanie przebiegów czasowych sygnału po filtrze zerującym i po wygładzeniu dla różnych parametrów a')%Nadanie nazwy wykresu

legend('Filtr Zerujacy', 'Wygładzenie dla a=0.7', 'Wygładzenie dla a=0.5', 'Wygładzenie dla a=0.9');%Ustawienie etykiet legendy

PoWyglFourier=fft(PoWygl2,CzestProb); %Transformata Fouriera sygnału po wygładzeniu

PoWyglGWM=abs(PoWyglFourier).^2./CzestProb; %Gęstość widmowa mocy po wygładzeniu

CzestNyq = CzestProb/2; %Obliczenie częstotliwości nyquista

figure; %Deklaracja nowego obszaru kreslenia

semilogy((0:CzestNyqZakl-1), PoZerGWM(1:CzestNyqZakl), 'Color', 'g'); %Rysowanie gęstości widmowej mocy sygnału EKG po filtracji filtrem zerującym

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

semilogy((0:CzestNyq-1), PoWyglGWM(1:CzestNyq), 'Color', 'r'); %Rysowanie gęstości widmowej mocy sygnału wygładzonego

grid on; %Włączenie siatki na wykresie

title('Porównanie gęstości widmowych mocy oryginalnego sygnału i po wygładzeniu')%Nadanie nazwy wykresu

legend('Filtr Zerujacy', 'Wygładzenie dla a=0.5'); %Ustawienie etykiet legendy

%Zadanie 6.

CzasKor=((-WymiarX/CzestProb+OkresProb):OkresProb:WymiarX/CzestProb-OkresProb)';%Obliczenie czasu autokorelacji

AutokorEKG = xcorr(SygnalEKG); %Autokorelacja oryginalnego sygnału EKG

AutokorZer = xcorr(SygnalEKGPoFiltrZer); %Autokorelacja sygnału EKG po filtracji filtrem zerującym

AutokorWygl = xcorr(PoWygl2); %Autokorelacja sygnału EKG po filtracji filtrem wygładzającym

figure; %Deklaracja nowego obszaru kreslenia

plot(CzasKor, AutokorEKG, 'g'); %Rysowanie wykresu autokorelacji oryginalnego sygnału

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(CzasKor, AutokorZer), 'r'; %Rysowanie wykresu autokorelacji sygnału po filtracji filtrem zerującym

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(CzasKor, AutokorWygl), 'b'; %Rysowanie wykresu autokorelacji sygnału po filtracji filtrem wygładzającym

grid on; %Włączenie siatki na wykresie

title('Porównanie przebiegów czasowych autokorelacji poszczególnych sygnałów')%Nadanie nazwy wykresu

legend('EKG', 'Filtr Zerujący', 'Wygładzenie'); %Ustawienie etykiet legendy

%Zadanie 7.

SygnalEKGBezOffsetu = SygnalEKG - mean(SygnalEKG); %Sygnał EKG bez offsetu

PoZerBezOffsetu = filter(FiltrZerujacyLicznik, FiltrZerujacyMianownik, SygnalEKGBezOffsetu);%Filtracja sygnału bez offsetu filtrem zerującym

PoWyglBezOffsetu = filter(FiltrWyglLicznik2, FiltrWyglMianownik2, SygnalEKGBezOffsetu);%Filtracja sygnału bez offsetu filtrem wygładzającym

PoprAutokorEKG = xcorr(SygnalEKGBezOffsetu); %Autokorelacja sygnału EKG bez offsetu

PoprAutokorZer = xcorr(PoZerBezOffsetu); %Autokorelacja sygnału EKG bez offsetu po filtracji filtrem zerującym

PoprAutokorWygl = xcorr(PoWyglBezOffsetu); %Autokorelacja sygnału EKG bez offsetu po filtracji filtrem wygładzającym

figure; %Deklaracja nowego obszaru kreslenia

plot(CzasKor, PoprAutokorEKG, 'g'); %Rysowanie wykresu autokorelacji sygnału bez offsetu

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(CzasKor, PoprAutokorZer), 'r'; %Rysowanie wykresu autokorelacji sygnału bez offsetu po filtracji filtrem zerującym

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(CzasKor, PoprAutokorWygl), 'b'; %Rysowanie wykresu autokorelacji sygnału bez offsetu po filtracji filtrem wygładzającym

grid on; %Włączenie siatki na wykresie

title('Porównanie przebiegów czasowych autokorelacji poszczególnych sygnałów bez offsetu')%Nadanie nazwy wykresu

legend('EKG', 'Filtr Zerujący', 'Wygładzenie'); %Ustawienie etykiet legendy

figure; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(CzasKor, PoprAutokorWygl), 'b'; %Rysowanie wykresu autokorelacji sygnału bez offsetu po filtracji filtrem wygładzającym

grid on; %Włączenie siatki na wykresie

axis([-2.5 , 2.5,-Inf ,Inf]); %Ustawienie limitów osi

title('Porównanie przebiegów czasowych autokorelacji poszczególnych sygnałów bez offsetu (Przybliżenie)')%Nadanie nazwy wykresu

legend('Wygładzenie'); %Ustawienie etykiet legendy

%Zadanie 5*.

K1 = 2; %Rząd średniej ruchomej 1

K2 = 10; %Rząd średniej ruchomej 2

K3 = 20; %Rząd średniej ruchomej 3

PoSredniejRuchomej1 = filter(ones(1,K1)/K1, 1, PoWygl2); %Filtracja sygnału przy pomocy średniej ruchomej 1

PoSredniejRuchomej2 = filter(ones(1,K2)/K2, 1, PoWygl2); %Filtracja sygnału przy pomocy średniej ruchomej 2

PoSredniejRuchomej3 = filter(ones(1,K3)/K3, 1, PoWygl2); %Filtracja sygnału przy pomocy średniej ruchomej 3

figure; %Deklaracja nowego obszaru kreslenia

plot(Czas, PoSredniejRuchomej1, 'g'); %Rysowanie przebiegu czasowego sygnału przefiltrowanego przy pomocy średniej ruchomej 1

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoSredniejRuchomej2, 'r'); %Rysowanie przebiegu czasowego sygnału przefiltrowanego przy pomocy średniej ruchomej 2

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoSredniejRuchomej3, 'b'); %Rysowanie przebiegu czasowego sygnału przefiltrowanego przy pomocy średniej ruchomej 3

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoWygl2, 'y'); %Rysowanie przebiegu czasowego sygnału wygładzonego

grid on; %Włączenie siatki na wykresie

axis([-Inf , 5,-Inf ,Inf]); %Ustawienie limitów osi

title('Porównanie przebiegów czasowych przed i po filtracji średnią ruchomą dla różnych wartości K')%Nadanie nazwy wykresu

legend('Średnia ruchoma dla K=2', 'Średnia ruchoma dla K=10', 'Średnia ruchoma dla K=20');%Ustawienie etykiet legendy

PoSredniejRuchomejFourier=fft(PoSredniejRuchomej2, CzestProb); %Transformata Fouriera sygnału po filtracji średnią ruchomą

PoSredniejRuchomejGWM=abs(PoSredniejRuchomejFourier).^2./CzestProb; %Obliczanie GWM

CzestNyquist = CzestProb/2; %Obliczenie częstotliwości nyquista

figure; %Deklaracja nowego obszaru kreslenia

semilogy((0:CzestNyquist-1), PoSredniejRuchomejGWM(1:CzestNyquist), 'Color', 'g');%Rysowanie gęstości widmowej mocy sygnału przefiltrowanego średnią ruchomą

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

semilogy((0:CzestNyquist-1), PoWyglGWM(1:CzestNyquist), 'Color', 'r'); %Rysowanie gęstości widmowej mocy sygnału wygładzonego

grid on; %Włączenie siatki na wykresie

title('Porównanie gęstości widmowych mocy przed i po filtracji średnią ruchomą')%Nadanie nazwy wykresu

legend('Średnia ruchoma dla K=10', 'Wygladzony'); %Ustawienie etykiet legendy

PoBandpass = filter(bandpass, PoWygl2); %Filtracja zaprojektowanym filtrem pasmowoprzepustowym

fvtool(bandpass); %Charakterystyki filtra pasmowoprzepustowego

figure; %Deklaracja nowego obszaru kreslenia

plot(Czas, PoBandpass, 'g'); %Rysowanie przebiegu czasowego sygnału przefiltrowanego przy zaprojektowanego w FDATool filtra pasmowoprzepustowego

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

plot(Czas, PoWygl2, 'r'); %Rysowanie przebiegu czasowego sygnału wygladzonego

grid on; %Włączenie siatki na wykresie

axis([8 , 16,-Inf ,Inf]); %Ustawienie limitów osi

title('Porównanie przebiegów czasowych przed i po filtracji filtrem pasmowoprzepustowym')%Nadanie nazwy wykresu

legend('Filtr Pasmowoprzepustowy', 'Wygładzony'); %Ustawienie etykiet legendy

PoBandpassFourier=fft(PoBandpass, CzestProb); %Transformata Fouriera sygnału po filtracji filtrem pasmowoprzepustowym

PoBandpassGWM=abs(PoBandpassFourier).^2./CzestProb; %Obliczanie GWM

figure; %Deklaracja nowego obszaru kreslenia

semilogy((0:CzestNyquist-1), PoBandpassGWM(1:CzestNyquist), 'Color', 'g'); %Rysowanie gęstości widmowej mocy sygnału przefiltrowanego filtrem pasmowoprzepustowym

hold on; %Dalsze rysowanie na tym samym obszerze kreślenia

semilogy((0:CzestNyquist-1), PoWyglGWM(1:CzestNyquist), 'Color', 'r'); %Rysowanie gęstości widmowej mocy sygnału wygładzonego

grid on; %Włączenie siatki na wykresie

title('Porównanie gęstości widmowych mocy przed i po filtracji filtrem pasmowoprzepustowym')%Nadanie nazwy wykresu

legend('Bandpass', 'Wygladzony'); %Ustawienie etykiet legendy

Wnioski:

Na podstawie przebiegu czasowego możemy stwierdzić, że zastosowanie średniej ruchomej za skutkowało wygładzeniem sygnału oraz zwiększeniem pików widocznych w oryginalnym sygnale.

Dokonałem porównania na jednym wykresie odpowiedzi czasowych z wykorzystaniem 3 różnych rzędów średniej ruchomej i piki są tym bardziej wzmacniane im większy jest rząd K.

Filtr pasmowo-przepustowy został zaprojektowany w środowisku FDATool tak aby przepuszczał częstotliwości z przedziału 0.5Hz do 40Hz gdyż to właśnie w tym paśmie znajduję się interesujący nas sygnał. Po jego zastosowaniu przebieg czasowy uległ przesunięciu lecz nie ma to większego wpływu na odczyt interesujących nas danych, ponieważ nieważne jest położenie pików w czasie, lecz ich wzajemna odległość na podstawie, której obliczamy puls. Filtr został utworzony poprzez wybranie poniższych opcji:

Fs = 250; % Sampling Frequency

Fstop1 = 0.2; % First Stopband Frequency

Fpass1 = 0.6; % First Passband Frequency

Fpass2 = 36; % Second Passband Frequency

Fstop2 = 40; % Second Stopband Frequency

Dstop1 = 0.001; % First Stopband Attenuation

Dpass = 0.057501127785; % Passband Ripple

Dstop2 = 0.001; % Second Stopband Attenuation

Wybrany przezemnie model to Filtr FIR, Window, Kaiser