Evaluating a raw ECG signal in MATLAB

Question

0 votes

Okay so I am loading a raw ECG signal. I know the duration of the signal is 690seconds. Using this information, I need to find the length of the signal, the sampling frequency, the energy, and the power. I need to create a filter, a bandpass. So here I created a lowpass and highpass... is there a way to combine them to get a bandpass? I'm lost on this. Next I used the function pwelch and try to plot my results. I need to take out frequencies above 10Hz (White Gaussian noise) and below 5Hz (Low background noise). Should the c in pwelch = (b-B)/(a-A) ?

This is the code I have right now, but I'm unsure if my power spectrum density graph is correct?.

    % Part 1
    % Loads raw ECG signal.
    load('pj3.mat');
    % Part 2
    % Duration = z = 690s. Length? Samp Freq? Energy? Power?
    L=length(data);                       %Length
    Fs=L/690 ;                            %Sampling Frequency           
    E=sum((data).^2);                     %Energy
    P=E/L;                                %Power
    cutoff1=10;                           %Low background noise
    cutoff2=5;                            %White Gaussian noise
    Wn1=cutoff1/(Fs/2);                   %Low Cutoff frequency
    Wn2=cutoff2/(Fs/2);                   %High Cutoff frequency
    [b,a]=butter(2, Wn1, 'low');          %2nd order lowpass filter
    [B,A]=butter(2, Wn2, 'high');         %2nd order highpass filter
    %Part 3
    [Pxx,f]=pwelch([c,[],[],[],(Fs)]);
    figure(1);
    plot(f,Pxx);                          %Plot power spectrum density
    title('Power Spectrum Denisty');
    xlabel('Frequency (Hz)')
    ylabel('Power(dB)')

    %Part 4
    % These are the values that are returned after part 2: A, B, a, b
    % A=[1,-1.955578240315036,0.956543676511203];
    % B=[0.978030479206560,-1.956060958413119,0.978030479206560];
    % a=[1,-1.911197067426073,0.914975834801434];
    % b=[9.446918438401619e-04,0.001889383687680,9.446918438401619e-04];
    n=0:pi/200:pi/2;
    %Plot HP mag/phase freq resp
    figure(2);
    freqz(B,A,n);                         
    title('Magnitude and Phase Frequency Response for High-Pass Filter');
    %Plot LP mag/phase freq resp
    figure(3);
    freqz(b,a,n);                         
    title('Magnitude and Phase Frequency Response for Low-Pass Filter');

    %Part 5
    %For High Pass Filter
    y1 = filtfilt(B,A,data);
    [pks1,locs1] = findpeaks(data);
    %For Low Pass Filter
    y2 = filtfilt(b,a,data);
    [pks2,locs2] = findpeaks(data);
    time=locs2-locs1;
    %To find accuracy:
    Anntest.time = time;
    Anntest.unit = ones( size(Anntest.time) );
    Anntruth = load_eaf('Anntruth.eaf');
    accuracy = eaf_compare(Anntruth, Anntest,'Window',0.15);

For part 5, am I calculating time right? I need to develop an approach to detect the peak of the R-wave. The final result should be a vector of times (in seconds) corresponding to the R-wave peak locations. Adding missing the R-wave manually is NOT allowed in this project. What functions will help me do this?