Signalerx4


ØVELSE 4
i 31610


2009-03-23

Rasmus Berg Palm s062241


Mads Andersen s062237


Jacob Bjerring Olesen s062248


Indholdsfortegnelse

Exercise

1

Et ECG signal har nogle karakteristika. Herunder en figur der benævner disse:

ECGintervals3.gif

En redegørelse for disse benævnelser kan findes her:
http://server.oersted.dtu.dk/31610/?exercises/exercise4.html

Vi er givet et ECG signal og skal forsøge at sige noget om det.
Nedenfor ses signalet som optaget ved leads (der er korrigeret for forstærkning)

atleads.png

Vi ligger her mærke til 3 ting:
Signalet har en periodisk komponent med en periode på lidt under 1s. Dette tænkes at være R til R-intervallet.
Signalet har en stor lav-frekvent komponent.
Signalet har en lille høj-frekvent komponent.
Herunder et zoom på en enkelt periode:

oneper.png

Der ses tydeligt P tak, QRS kompleks samt T tak.
Da vi skal forsøge at sige noget kvalitativt om disse vil vi forsøge at øge signal-støj forholdet ved at fjerne den lav og høj-frekvente komponent.
Vi konstruerer nu et ideelt højpasfilter i frekvensdomænet og plotter det sammen med spektret af signalet:

transferfunction.png

Vi vil altså fjerne alt under 0.5Hz med dette filter.
For at sikre os at dette filter har et reel impulsrespons må vi kigge på det i tidsdomænet. Ved en ifft fås:

impulseresponse.png

Det ses at impulsresponset er reelt. Dette var forventet da spektret var en lige funktion af f (2.12).
Filtreringen foregår i frekvensdomænet ved at gange overføringsfunktionen H med spektret for signalet og derefter ifft'e det filterede signal tilbage til tidsdomænet:

highfiltered.png

Det ses at den lavfrekvente komponent er forsvundet. Lad os se hvordan denne så ud. Dette gøres ved at plotte forskellen mellem det filtrerede og det ufiltrerede signal:

lowfreqpart.png

Dette er næsten for pænt til at være biologisk. Vi mistænker overlagt forurening af signalet.
Følgende matlab kode blev brugt til ovenstående:

% Ex1
clear all; close all; clc;
load ecg

% sampling frequency:
fs = 500; %Hz

% time axis:
t = 0:1/fs: (length(ecg)-1)/fs;

% the signal before amplification:
ecg = ecg/500;

% the signal measured at the leads:
figure;
plot(t, ecg)
xlabel('time (s)')
ylabel('voltage (V)')
title('The ECG signal as measured at the leads')

% The spectrum of the signal:
M = length(ecg);
ECG = fftshift(fft(ecg)); %vi laver ikke zeropadding, for sådan ser signalet ikke ud...(men man kan sige at vi antager det er periodisk...)
F = [0:M-1]/M*fs-fs/2;
ECGabs = abs(ECG);

% The filter transfer function:
H = zeros(size(ECG));
f0 = 0.5;
H(F>=f0)=1;
H(F<=-f0)=1;

% Plot of transfer function and the spectrum:
figure;
plot(F,ECGabs, F, H);
legend('|ECG spectrum|', '|H(F)|')
title('Plot of transfer function and the spectrum of the unfiltered signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude (V)')
axis([-10 10 0 1.1])

% The impulse response of the filter:
h = ifft(ifftshift(H));
% dt = 1/fs;
% n = 0:length(h)-1;
% ha = -(sin(2.*pi.*f0.*n.*dt)-sin(pi.*n))./(pi.*n);
% ha(1) = -2*f0*dt+1;
figure;
plot(t,h)
title('The impulse response of |H(F)|');
xlabel('Time (s)');
ylabel('Amplitude (V)');
% impulsresponset er reelt, da overføringsfunktionen er lige!

% The filtering is done in the frequency spectrum:
ECGfilt = ECG.*H;

% Inverse fourier transforming:
ecgfilt = ifft(ifftshift(ECGfilt));

figure;
plot(t,ecgfilt)
title('The hihgpass filtered signal');
xlabel('Time (s)');
ylabel('Amplitude (V)');

figure;
plot(t,ecg-ecgfilt)
title({'The difference between the hihgpass-filtered and the un-filtered signal' '(The low frequency component of the unfiltered signal)'});
xlabel('Time (s)');
ylabel('Amplitude (V)');

2

Nu til den højfrekvente støj. Lysnettet giver som regel en vis mængde støj ved 50Hz. Denne vil vi nu forsøge at filtrere fra ved at benytte et notch filter omkring 50Hz. Da vi forventer at støjen er fra virkeligheden og altså ikke ligger perfekt i 50Hz benytter vi et notch filter fra 49Hz til 51Hz. Vi konstruerer et ideelt filter i frekvensdomænet:

notchfilter.png

Og filtrerer ligesom før i frekvensdomænet:

highnnotchfiltered.png

Der ses en mindre forbedring af signal-støj forholdet, mest synligt i TQ intervallet hvor hjertet ikke har den store elektriske aktivitet.
Følgende matlab kode blev brugt til ovenstående:

% Ex2
% The filter:
Hnotch = ones(size(ECGfilt));
Hnotch(F<51 & F>49)= 0;
Hnotch(F>-51 & F<-49)= 0;

figure;
plot(F, Hnotch)
title('The notchfilter transfer function')
xlabel('Frequency (Hz)')
ylabel('Amplitude (V)')
ylim([0 1.1])

% The filtering:
ECGfilt2 = ECGfilt.*Hnotch;

% The filtered signal in the time domain:
ecgfilt2 = ifft(ifftshift(ECGfilt2));

figure;
plot(t,ecgfilt2)
title('ecg signal filtered with notch, and highpass filters')
xlabel('Time (s)');
ylabel('Amplitude (V)');

3

Selv efter at signalet er filtreret to gange, først for muskel støj og dernæst for det elektriske netværk, vil det stadig være påvirket af andet støj. Forholdet mellem det brugbare signal og støjen kan dog varieres alt efter hvilken båndbrede filteret har. Mængden af støj er proportionalt med bredden af båndet, det vil sige at for små båndbredder er støjen tilsvarende lille, dog er der risiko for at fjerne brugbar information fra signalet, hvilket ikke er hensigtsmæssig.
Det bedste signal-støj forhold er fundet ved en cut-off frekvens på omkring 27 Hz. Signalet for dette kan ses i nedenstående figur.
opg3-27.png

Plottes der for en cut-off frekvens på 15, ses det at signalet bliver betydelig svagere, samtidig med at det er sværere at skelne imellem de forskellige intervaller af ECG'et.
opg3-15.png

Følgende figur viser de forskellige intervaller af ECG'et
opg3-27i.png

% Ex3, the cut-off frequency comprimise is found to 27,
H=zeros(size(ECGfilt2));
H(F<27 & F>-27) = 1;
ECGfilt_27 = ECGfilt2.*H;
ecgfilt_27 = ifft(ifftshift(ECGfilt_27));

t = [0:length(ecgfilt_27)-1]/fs;

figure;
plot(t,ecgfilt_27);
xlabel('Time (s)')
ylabel('Voltage (V)')
title('cutoff frequency at 27')
ylim([-5e-4 12e-4]);

H=zeros(size(ECGfilt2));
H(F<15 & F>-15) = 1;
ECGfilt_15 = ECGfilt2.*H;
ecgfilt_15 = ifft(ifftshift(ECGfilt_15));

t = [0:length(ecgfilt_15)-1]/fs;

figure;
plot(t,ecgfilt_15);
xlabel('Time (s)')
ylabel('Voltage (V)')
title('cutoff frequency at 15')
ylim([-5e-4 12e-4]);

4

ecg signalet er tilnærmelsesvist periodisk, idet hjetet jo netop trækker sig sammen periodisk, med en puls der svarer til 1/perioden. Dette udnyttede vi også da vi fouriertransformede tidligere, idet fft således gav en fornuftig beskrivelse i frekvensdomænet uden zeropadding af signalet. Perioden, og dermed pulsen kan findes ved autokorrelation af signalet - de tidsforskydninger hvor signalet korrelerer bedst med sig selv er et heltal gange perioden. Da signalet ikke er perfekt periodisk, vil den største top i autokorrelationssekvensen naturligvis være i tidsforskydning 0. Perioden findes derefter ved den næste 'største' top, eller som afstanden mellem to nærmeste 'største' toppe. Da signalet har endelig længde, vil toppene som markerer perioden dog aftage i størrelse som tidsforskydningen øges, da der 'zeropaddes'. Den næststørste top i autokorrelationssekvensen vil således være den første gang at signalet gentager sig selv. Den næststørste top er altså lig den næste 'største' top, efter den i tidsforskydning 0. Hvis man bruger MATLABs max funktion og søger på et område der ikke inkluderer toppen i tidsforskydning 0, men alle de andre toppe, er man således sikker på at det er den top i tidsforskydning=1*perioden som findes. Funktionen heartrate.m er konstrueret til at finde perioden og pulsen af et hjerte ud fra et ecg signal. For at mindske beregningsantalellet søges kun efter pulsen i området 30bpm til 220bpm. Den nedre grænse ville sagtens kunne flyttes til en lavere værdi uden at det gav et andet resultat (periode=1/puls, så nedre grænse puls er øvre grænse periode), men den øvre grænse kan ikke gøres meget højere, da man kommer for tæt på toppen i tidsforskydning 0. Det burde dog heller ikke være nødvendigt at søge i et interval med en højere max-puls, da den fysiologiske max-puls for mennesker ca. er 220-alder. Nedenfor ses heartrate.m.

function [hr, per, ecgcorr, lags] = heartrate(ecgsig,fs, doplot)
% Heartrate finds the heartrate from and ecg signal. If there is more than
% 2 input arguments, heartrate plots  the autocorrelation function of the 
% ecg signal.

% The autocorrelation of the signal:
[ecgcorr,lags] = xcorr(ecgsig,'coeff');

% lags in seconds:
lags = lags/fs;

% Heart rate limits:
lowlimr  = 30; %bpm
highlimr = 220; %bpm

% Normal heart period limits (search interval):
highlimp  = 1/lowlimr*60; %s 
lowlimp = 1/highlimr*60; %s

% Making the search interval window:
ecgcorrW = ecgcorr(lags>lowlimp & lags<highlimp);
lagsW = lags(lags>lowlimp & lags<highlimp);

% Finding the period of this heart in vector indices:
[maks,peri] = max(ecgcorrW);

% The corresponding lag in seconds (heart period in seconds):
per = lagsW(peri);%s

% The heart rate:
hr = 1/(per/60);%bpm

% plot of autocorrelation function:
if nargin > 2
figure;
plot(lags,ecgcorr)
xlabel('lags (s)')
ylabel('Normalized autocorrelation sequence')
title('Autocorrelation of ecg signal')
end

Heartrate.m er brugt på det filtrede, hhv. rå ecg-signal, hvilket gav følgende værdier og autokorrelationssekvenser:
Heart period from ecgfilt2: 0.912s
Heart rate from ecgfilt2: 65.7895bpm

Heart period from unprocessed ecg signal: 0.274s
Heart rate from unprocessed ecg signal: 218.9781bpm
autofilt.png
autoufilt.png

Perioden fundet ud fra det filtrede signal passer godt med afstanden mellem eksempelvis to R-takker i ecg-signalet, og må derfor antages korrekt. Autokorrelationssekvensen ser ud som spået, nemlig med 'største' toppe med lige store mellemrum, og aftagende i størrelse som tidsforskydning øges. Det ufiltrede signal har ikke samme periode, pga. de forurende lavfrenkvenser, og autokorrelation bliver derfor helt anderledes. Det er umuligt umiddelbart at aflæse den korrekte periode ud fra det ufiltrede signals autokorrelationssekvens, og pulsen der findes er også helt forkert.
Ud over heartrate.m er følgende MATLAB-kode konstrueret til denne opgave:

% Ex4
% Finding heartrate from the processed signal ecgfilt2:
[hrfilt2, perfilt2, corrfilt2, lagsfilt2] = heartrate(ecgfilt2,fs);
display(['Heart period from ecgfilt2: ' num2str(perfilt2) 's'])
display(['Heart rate from ecgfilt2: ' num2str(hrfilt2) 'bpm'])

%figure of autocorrelation function:
figure;
plot(lagsfilt2, corrfilt2)
xlabel('lags (s)')
ylabel('Normalized autocorrelation sequence')
title('Autocorrelation of ecgfilt2')
axis([-8.4 8.4 -0.3 1.1])

% Finding heartrate from the unprocessed signal ecg:
[hr_un, per_un, corr_un, lags_un] = heartrate(ecg,fs);
display(['Heart period from unprocessed ecg signal: ' num2str(per_un) 's'])
display(['Heart rate from unprocessed ecg signal: ' num2str(hr_un) 'bpm'])

%figure of autocorrelation function:
figure;
plot(lags_un, corr_un)
xlabel('lags (s)')
ylabel('Normalized autocorrelation sequence')
title('Autocorrelation of ecgfilt2')
axis([-8.4 8.4 -0.5 1.1])

edit sections

Beskeder
Click me to edit !
Drag me !
Click me to edit !
sigx4
Here is a place for your title Click me to edit ! false
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License