How do you repeat a QRS signal?

I have a code that currently plots the magnitude and phase components of an ECG signal. Currently I have it plot only one cycle of a wave in Figure 1, but I need it to plot more than one so I can calculate the time duration between the two consecutive R peaks. I also want to see if I can determine the bandwidth of the signal using 'bw = powerbw(x)' command. Finally, is there a way to lanel the frequencies with the largest magnitude?

clear all; close all;
load ECG1; 
fs = 125; 
N = length(ecg);
t = (1:N)/fs;
plot(t,ecg); 
xlabel('Time (sec)');
ylabel('ECG');
figure;
f = (1:N)*fs/N;
N2 = round(N/2);
X = fft(ecg);
Mag = abs(X);
Phase = (unwrap(angle(X)))*360/(2*pi);
subplot(2,1,1);
plot(f(1:N2-1),Mag(2:N2));
xlabel('Frequency (Hz)');
ylabel('|ECG(f)|');

subplot(2,1,2);
plot(f(1:N2-1),Phase(2:N2));
xlabel('Frequency (Hz)');
ylabel('Phase ECG (deg)');

Answer 1

[William Rose]

@collegestudent,

If you want to use frequency analysis to examine the heart rate frequencies, you will need a lot of beats. 300 seconds is the minimum re ommended duration for analysis of heart rate variability.

The ECG has fast components, associated with the QRS complex, and slow components associated ith heart rate variability.

If you want to use frequency analysis to examine the heart rate and its variability, you might want to check out this paper: Spectral analysis of point processes, Rose & Knight, Med.Sci.SportsExercise 2011.

Can you upload your file, ECG1(.mat)?

[William Rose]

@collegestudent,

[edit: fix typos]

I suspect from your question and your code that you are relatively new to frequency analysis. But you are doing it - you are ahead of where I was, when I was in college.

Fourier analysis of a finite length signal creates the spectrum that will reproduce the signal if the sample is one cycle of a periodic signal. Therefore, if the beginning and the end have different values or slopes, the spectrum will include a bunch of high frequency power to handle the wrap-around discontinuity at the end/start. Which is why people use tapering windows, such as Hamming, Hann, etc. Look it up.

Your code plots all the frequencies, including the "top half" of the spectrum (frequencies from fs/2 to fs). This is unnecessary, since the top half of the spectrum is just a reflected copy of the lower half spectrum.

Yes it is possible to label a point on the spectrum, such as the max value. Use

[MagMax,idx]=max(Mag);

to find the max value of the amplitude spectrum. Then the associated frequency is

fmax=f(idx-1);  

where the "-1" is because your frequency vector is offset by 1 from the Mag vector. Then use text() to place the desired label on the plot - see help for text(), which includes good examples. For example, you could do

textStr=sprintf('Max power at %.2f Hz',fmax);   %create string
text(fmax,MagMax,textStr);                      %place string on plot

Good luck.

[collegestudent]

@William Rose I have added ECG1.mat

[collegestudent]

So if i want to label the three largests magnitudes would I have to do

[MagMax,idx]=max(Mag);
[MagMax2,idx]=max(Mag);
[MagMax3,idx]=max(Mag);

[William Rose]

@collegestudent,

[magMax,idx]=maxk(Mag,3);

returns the three largest values and their indices, as vectors.

[William Rose]

@collegestudent,

I have looked at the signal. Remove the mean value before computing the spectrum to avoid non-useful power at and hear f=0. I see now that the file includes a 1-beat-long record and a two-beat-long record.

[collegestudent]

When I used [magMax,idx]=maxk(Mag,3); I got an error saying that MAX with two matrices to compare and two output arguments is not supported.

[William Rose]

@collegestudent,

ecg1=load('ECG1');
fs=ecg1.fs;  
ecg=ecg1.ecg-mean(ecg1.ecg);    % remove mean value
N = length(ecg);
N2 = round(N/2);
t = (1:N)/fs;                   % time vector
f = (0:N-1)*fs/N;               % frequency vector
Mag = abs(fft(ecg));            % magnitude spectrum
%% Plot results
plot(t,ecg,'-r.'); 
xlabel('Time (sec)'); ylabel('ECG'); grid on

figure;
plot(f(1:N2),Mag(1:N2),'-r.');
xlabel('Frequency (Hz)'); ylabel('|ECG(f)|'); grid on
[magMax,idx]=maxk(Mag(1:N2),3);
for i=1:3
    textStr=sprintf(' %.2f Hz',f(idx(i))); % create string
    text(f(idx(i)),magMax(i),textStr);    % place string on plot
end

The 3 highest-power frequencies are labelled.

Try it.

[collegestudent]

Thank you for this. Is there a way to determine bandwidth of signal and average heart rate using MATLAB?

[William Rose]

@collegestudent,

Yes you can analyze heart rate data in the frequency domain with Matlab.

To get average heart rate, simply count the number of R waves in a given duration, such as 1 minute. Search Matlab Answers and the File Exchange for how to find the R waves. Learn how to use findpeaks() by reading the help and trying the examples.

The term "bandwidth" is not generally applied to heart rate because it is not useful.

The EKG includes fast and slow components, which convery different types of information. The fast components create the wave shape details: S-T segment height; duration of intervals such as P-R, Q-T, etc. These details are essential for diagnosing many types of condtions, when done with a 12 lead EKG. The frequency domain analysis of the raw EKG has not yielded important clinical insights or tools. The analysis is best done on the time domain signal.

The slow components of the EKG (slower than the HR frequency) pertain to how HR changes with time. This is referred to as HR variability. For this kind of study, you typically collect data for 5 minutes or more. You determine the beat times (the R wave times). There are many research articles about HR variability. See here for an introduction. See here for how to compute the spectrum.