1

signal processing is not my domain of expertise , but I need to make use of it, to do an HRV analysis. These are steps I do starting from my raw ECG signal.

ECG signal -> RR peak detection -> find RR intervals -> cubic spline interpolation -> IIR high pass filter -> FFT

The peaks I get is a mirrored image, I guess its because its a pure signal from a digital patient simulator like all at 60 bpm or all at 120 bpm etc and follows the rule of Hermetian symmetry.

However my FFT output at pure 60 bpm should have indicated a peak at 1 Hz and a 120bpm at 2 Hz, is what I understood. But my all peaks start from 0.

Please help

PS : Implementing this all in an android application.

PS : This is what my output looks on a 60 bpm pure signal.

Three graphs in order : RR interval, RR after interpolation, FFT

PS : Added log data

FFT output : http://www.yourfilelink.com/get.php?fid=1337132

Diwesh Saxena
  • 11
  • 1
  • 3

1 Answers1

1

The FFT output will be bins evenly spaced from DC (bin 0) to 1 bin less than your sampling frequency. Thus your sampling frequency is at N for bins 0 to N-1, where N is the number of samples in your time domain sequence.

Within the first Nyquist zone of $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate, the bins from 0 to N/2-1 are the positive frequency components, while the bins from N/2 to N-1 are the negative components. Specifically bin N-1 is -1, bin N-2 is -2, etc so you could view the domain as a circle from 0 to N-1 and repeating again.

The reason your peaks are a mirror image is because your time domain signal is real, and for real signals the positive and negative frequencies have a conjugate symmetric relationship (same magnitude, opposite phase).

The reason you see a large response in bin 0 is because your time domain signal has a non-zero mean. Bin 0, being the DC term, is the mean of your signal (scaled by a factor of N). You likely don't care about DC, and would therefore benefit by removing the mean of your signal before taking the FFT.

The FFT is an algorithm that implements the DFT efficiently, which is as the name implies a discrete version of the Fourier Transform. I think it is helpful to see that the Fourier Transform is simply a correlation as we sweep frequency. (Where correlation is complex conjugate multiply and integrate, or equivalently discrete as complex conjugate multiply and accumulate). Thus for the DFT, each bin is a correlation of a single rotating phasor (your traditional sine wave is made of two equal and opposite rotating phasors as shown by Euler's Identity: $\frac{1}{T}(e^{j\omega t}+e^{-j\omega t}$). So for the DFT, given as: $$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j(k\omega_0) n}$$

So we see that for the first bin (k=0), we just have N times the mean of x. If x has a DC value, we would get a strong correlation for that calculation. Also observe when k =1, which corresponds to the next bin, the exponential is spinnning at its lowest possible rate $\omega_0$ corresponding to one cycle over your time domain sequence. So in that calculation we are computing the correlation of the lowest possible frequency that is spinning in the opposite direction ($e^{+j\omega t}$): such a signal in our x(t), after this computation would have the same result as bin 0: we would get N times the mean after having "derotated" the signal that was present there. Thus every bin goes through the same "averaging filter" which is a Sinc function in frequency (not a great filter), so our output at any given bin will be sensitive to some degree of frequencies everywhere in our spectrum (and as we see with the Sinc function there will be discrete frequencies that are at null positions and therefore will be locations where there is no "leakage" into other bins).

This should give you enough insight into some important fundamental points to dig further into the DFT/FFT. In addition here are some other related posts I think will help you:

Frequency measurement from acoustic emission sensor data(VDC RMS)

What proportion of a padded FFT should be actual values

Signal Processing/FFT gives very high magnitudes for low frequencies

You will want to read up a little more about windowing too; the links above and other posts (by search) include more information on that. Ultimately fred harris' paper has been a great resource for me on that topic.

Dan Boschen
  • 50,942
  • 2
  • 57
  • 135
  • Thanks for the useful resources and great explanation. So when you say I should essentially remove the DC components before applying the FFT, are you indirectly saying applying a high pass filter before FFT ? Does it do the same thing ? Also I added my output please see – Diwesh Saxena Mar 24 '17 at 00:04
  • Not really, what I suggest is actually a bad way to do a high pass filter but is very simple: the first bin is this mean value so we are essentially removing that; but generally not the way I would recommend filtering. (There are other posts that detail that but is essentially the "frequency sampling" method of filter design – Dan Boschen Mar 24 '17 at 00:10
  • Hmm, well since my knowledge is limited, I did try to apply an IIR second order high pass filter and then passed the output of that to FFT. As you can see from my graph the result of a 60 bpm data attached, doesnt seem right.Looking at the image can you think of probably where did I go wrong ? – Diwesh Saxena Mar 24 '17 at 00:26
  • Try my method if you didn't yet as although not ideal will be a little more fool proof – Dan Boschen Mar 24 '17 at 00:49
  • I did , removing the first bin and then plotting. bt as you can see from the graph I attached, even if i remove the first half, the graph still plots from max peak at 0, which is why I am so confused :( – Diwesh Saxena Mar 24 '17 at 00:54
  • I see.... I can't look at it until later unfortunately – Dan Boschen Mar 24 '17 at 01:07
  • I appreciate. I can post my IIR filter code too.I heard the coefficients chosen matter the most. Do you think posting the filter code will be helpful ? – Diwesh Saxena Mar 24 '17 at 01:10
  • Not sure yet - peaking at your plots briefly in between other things- is your horiz axis freq? A log (dB) magnitude would be more useful – Dan Boschen Mar 24 '17 at 01:20
  • yes it is. Yeah my app collects all data from RR to FFT. Ill try posting the text files – Diwesh Saxena Mar 24 '17 at 01:26
  • @DiweshSaxena - I looked at your data in log scale, and do see how you elminated the mean. I did not see from anything you wrote what the sampling rate is but see that you have 2^15 samples. If you want to be able to see 1 and 2 Hz signals if your FFT your data length in time will need to be >1sec, and really likely greater than 4 seconds to distinguish frequencies at 1Hz from DC. (as frequency resolution is 1/T). If your sampling rate for the data you provided as less than 8KHz we may be ok. The next thing to consider is windowing the data before taking the FFT – Dan Boschen Mar 24 '17 at 05:20
  • so the sampling rate of ecg device is 500 Hz. and I think there was this sampling rate variable I was applying to during interpolation which wsa 1000d. – Diwesh Saxena Mar 24 '17 at 05:42