How can I automatically identify and remove a fundamental frequency and all of its harmonics?

Question

I have some microscopy data that is contaminated by a heartbeat artifact that I'd like to remove. The data consists of a large time series of images captured at ~60Hz.

Here's a small example clip in GIF form:

I have taken the average pixel intensity over time, and computed the periodogram using Welch's method:

As you can see there is a sharp peak at ~1.8Hz which is likely to correspond to the heart rate (~108 beats/min). There are also a bunch of harmonic peaks at integer multiples of 1.8Hz. The exact heart rate is likely to vary from dataset to dataset, but I can specify a biologically plausible range as shown by the shaded area on the periodogram.

What I'd like to be able to do is:

1. Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics
2. Filter the data so as to remove the fundamental and all harmonics.

At the moment I can solve point 1 very crudely by finding the largest peak in the periodogram, then multiplying it by $1, 2, ..., N$ where $N$ is the estimated number of harmonic peaks, but I'm sure that there must be a better method than this hack.

Regarding point 2, I came across this question which mentions using a comb filter to remove a fundamental and all of its harmonics. Is this the best method to use? One important consideration is that I will have to apply the filter to each pixel timeseries in a large array, so a computationally efficient method would be highly desirable.

Example data

• I've uploaded a 1D timeseries of example mean pixel values in .csv format here.
• I've also added a sample of 1000 frames (spatially decimated 2:1) in a compressed .mat file here

Answer 1

Your method isn't that bad for a first attempt.

However, the following method tends to work better:

1. Search for local maxima
2. Check for maxima that are close (2 or 3 bin spacing) and merge them
3. Create a few hypotheses about the fundamental frequency. You currently assume the highest peak is the fundamental frequency, which is one hypothesis. You should also check for the possibility that the highest peak is the first harmonic, i.e. that there's a smaller peak at half the frequency of the main peak. You may also have other cases to consider, using your knowledge of the problem at hand (Interference? Irregular heartbeat?).
4. Assuming each of those hypothesis, find the ground frequency by fitting parabola to each harmonic peak. Each peak will produce a slightly different estimate due to noise, but these errors are uncorrelated and average out. One of the hypotheses will lead to much better fits, pick the ground frequency predicted by this hypothesis.
5. Using the ground frequency you found in step 4 as a given, refit the parabola around each peak to estimate the height of the peak. Note that the peaks probably fall between two bins.
6. You now have the location and strength of the fundamental and its harmonics, but not the phase. It's probably easiest to find the phase of the fundamental, subtract that, find the phase of the first harmonic, etc.

The core reason that this works better is step 4. Any wrong hypothesis for a ground frequency will fail horribly as you're trying to fit peaks around the location of predicted harmonics. Say you've got a peak at 2 Hz. This could be the ground frequency or a first harmonic. When you test the "first harmonic hypothesis", i.e. whether the ground frequency is actually 1 Hz, you fit parabola's to the data around 1,2,3,4,5 ... Hz. If that hypothesis is wrong, you get garbage around 1,3,5 Hz. If it's right, you might find peaks near 1.1 Hz, 2.2, 3.3, 4.4 and 5.6 - which suggests that the actual ground frequency is 1.12 Hz.

Answer 2

You're looking for Iterative Spectral Subtraction. Here is some general info from a Content Analysis book by Alexander Lerch.

Answer 3

I would recommend an auto-correlation on your periodogram. You can construct notch or negative peaking filters based on multiples of the fundamental frequency that the correlation yields.

This code helped me create a good auto-correlation plot (cepstrum is another good method when the harmonics have more power than the fundamental) http://note.sonots.com/SciSoftware/Pitch.html

Answer 4

1. Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics

You can take a local average across the DFT to find the heartbeats. If a certain point in this group is greater than threshold and the greatest value in some range around it then it is a heartbeat or harmonic.

2. Filter the data so as to remove the fundamental and all harmonics.

You could just use the previously located indices of heartbeats and stitch across them linearly. If that's not too crude for you.

What is the purpose of this part of the operation? Is this a medical analysis that requires the absence of the heartbeat to look for other (perhaps) hidden data? Or is this par just for aesthetics?

I imagine a comb filter will be somewhat uncontrollable. It will require careful adjustment of the feedback parameter. It could be viable.

EDIT: You require a pitch detection algorithm? I wrote one of these a few years ago by ignoring the (close to) DC terms and finding the first local maximum. Then using the bins on either side I could use a quadratic interpolation and find the location of the maximum frequency more accurately than the resolution allowed by the frequencies of the bins themselves.

Would it be more accurate to utilise harmonics? It would probably just add an opportunity for errors. Though you can try getting a value for the fundamental with this method and finding the maximum at double, triple etc and then using a similar method as before:

1. find local maximum 2. quadratic interpolate to find inter-bin value for frequency.

If you veto the quadratic step and just take the local max you will not get an accurate frequency and moving up the harmonics would help.