Power Spectrum of an ECG signal

We are going to apply the concepts of the previous section to the calculation of the power spectrum of a real electrocardiographic signal. We begin by loading the signal

 
-->stacksize(20000000)
 
-->load signals;
 
-->ecg = signals(12);

We can see that the sampling rate was 100 Hz

-->duration
 duration  =
 
    10.  
 
-->nsamples(12)
 ans  =
 
    1000.  
 
-->duration/nsamples(12)
 ans  =
 
    0.01

From an EDF file we can extract information to scale the signal with the next code

-->physmin(12)
 ans  =
 
  - 1920.2  
 
-->physmax(12)
 ans  =
 
    1926.  
 
-->digmin(12) 
 ans  =
 
  - 2048.  
 
-->digmax(12) 
 ans  =
 
    2047.  

-->physdim(12)
 ans  =
 
 uV

The result means that a value of -2048 corresponds to -1926 uV and a value of 2047 to 1926 uV. From these values, the amplitude of the signal can be corrected. But now we are going to calculate the power spectrum from the original values.

We are going to establish an epoch of one second corresponding to 100 values. We are going to build a vector containing the frequency scale

-->inct = 0.01
 inct  =
 
    0.01  
 
-->t = [1:100]*inct;
 
 
-->incf = 1/(100*inct)
 incf  =
 
    1.  
  
-->fre = [0:99]*incf;

Now we are going to consider three epochs and calculate the spectrum of each one of them

-->ndx = 1:100;                             // first epoch
 
-->x1 = ecg(ndx); 
 
-->xsetech([0,0,1/2,1/3])
 
-->plot2d(t,x1)
  
-->xsetech([1/2,0,1/2,1/3])
 
-->ffx1 = fft(x1,-1);
 
-->plot2d(fre,abs(ffx1)^2 /100);
 
-->ndx = 101:200;                           // second epoch
 
-->x2 = ecg(ndx);
 
-->xsetech([0,1/3,1/2,1/3]) 
 
-->plot2d(t,x2)
  
-->xsetech([1/2,1/3,1/2,1/3])
 
-->ffx2 = fft(x2,-1);
 
-->plot2d(fre,abs(ffx2)^2 /100); 
 
-->ndx = 201:300;                           // third epoch
 
-->x3 = ecg(ndx);
 
-->xsetech([0,2/3,1/2,1/3])
 
-->plot2d(t,x3)
 
-->xsetech([1/2,2/3,1/2,1/3])
 
-->ffx3 = fft(x3,-1);
 
-->plot2d(fre,abs(ffx3)^2 /100);

This result shows that even although the signals are very similar, the power spectrum varies considerably. Notice that frequencies above 50 Hz are above Nyquist frequency and that the spectra are symmetrical

Figure 12.3: Three epochs of `ecg' and its power spectra.

Now we are going to calculate the spectrum by using `pspect' with the next code

-->xbasc()
 
-->xsetech([0,0,1,1])
 
-->Pxx = pspect(50,100,'re',ecg); 
 
-->plot2d(fre,Pxx);

`Pspect' calculates the power spectrum of the whole signal after dividing it into epochs of 100 values. Notice that this spectrum is smoother than the individual spectra of the previous figure.

Figure 12.4: Power spectrum of `ecg' calculated using `pspect'

At this moment, we could ask whether the use of a longer epoch would have produced the same result. The next figure clearly shows that it is not the case. With the next code we are going to calculate the spectrum with epochs progressively longer

-->xbasc()
 
-->n = 200;                          // epoch of 200 values
 
-->fre  = [0:n-1] * (1/(n*inct));
 
-->ffx1 = fft (ecg(1:n),-1);
 
-->xsetech([0,0,1,1/3])
 
-->plot2d(fre, (abs(ffx1)^2)/n);
 
-->xtitle('n = 200')
 
 -->n = 500;                        // epoch of 500 values
 
-->fre  = [0:n-1] * (1/(n*inct));
 
-->ffx2 = fft (ecg(1:n),-1);
 
-->xsetech([0,1/3,1,1/3])
 
-->plot2d(fre, (abs(ffx2)^2)/n);
 
-->xtitle('n = 500')
  
-->n = 1000;                        // epoch of 1000 values
 
-->fre  = [0:n-1] * (1/(n*inct));
 
-->ffx3 = fft (ecg(1:n),-1);
 
-->xsetech([0,2/3,1,1/3])
 
-->plot2d(fre, (abs(ffx3)^2)/n);
 
-->xtitle('n = 1000')

Notice that, by increasing the length of the epoch, we increase the definition in frequency but we do not eliminate the variability among frequencies (increment of frequency becomes smaller). We need to average the epochs to get a clean, smooth power spectrum (alternatively we could average the values inside the spectrum).

Figure 12.5: Power spectrum of an ECG signal with epochs of different length. Notice that increasing the length of the window increases frequency resolution but do not smoothes the spectrum.