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

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
• Assuming that we're talking about a sequence of microscopy images in form of an animated picture, is there any chance to see that material? There may be a number of simpler and more accurate methods based on that signal if there is an above average spatial coherence in the artefacts. Reducing to the average pixel values will then throw away valuable information. – Jazzmaniac Feb 16 '15 at 14:06
• @Jazzmaniac I've added a GIF showing a small clip from an example dataset - there is definitely some local spatial coherence. – ali_m Feb 16 '15 at 15:43
• @Jazzmaniac I've also uploaded a sample frame sequence - I can provide more example data if it would be helpful, but the raw stacks themselves are much too big to upload. – ali_m Feb 16 '15 at 17:11
• You might be able to use a (musical) pitch detection/estimation method to estimate the fundamental frequency. – hotpaw2 Feb 17 '15 at 2:23
• @Jazzmaniac It's true that the heartbeat artifact is subtle to the human eye, but I don't really care that much about how the frames look - I'm interested in much more subtle local changes in intensity for which the heartbeat is a problem. The high-frequency oscillation that gives rise to the banding you see in single frames is a different issue with the recording equipment, which I planned to deal with separately. – ali_m Feb 18 '15 at 3:30

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.

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

• Thanks for the references. My first impression is that the speech enhancement methods are probably overkill given that the artifact I'm trying to remove is periodic and has a single fundamental frequency. I think I want to be looking more along the lines of pitch detection. – ali_m Feb 18 '15 at 3:57

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

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.

• The first part sounds more or less like the approach I suggested in my question, except your proposed method would not take into account the fact that harmonic frequencies should be regularly spaced. I don't really understand what you mean by "use the previously located indices of heartbeats and stitch across them linearly" - could you elaborate? The filtering step is necessary in order to analyse subtle local changes in intensity over time, rather than for aesthetic purposes. I've uploaded a slightly smaller .gif ("small" is a relative term here!). – ali_m Feb 18 '15 at 3:39
• 1. Comb filters are by definition: regularly spaced. 2. The issue is that adaptively filtering out the heartbeats with a comb filter of the fundamental freq is going to be hand-wavy and the result will be hand-wavy. So I don't think you could really rely on the data at these points too much. 3. Can you use different parts of the image and subtract them to destructively cancel out the heart beat, but maintain the data of interest? – Andrew Gallasch Feb 18 '15 at 5:46
• 1. My point was that the positions of the harmonic peaks should contain information that could be used to more accurately identify the fundamental frequency. 2. What exactly do you mean by "hand-wavy"? 3. That's something I hadn't considered. Nothing obvious occurs to me, but I'll think some more on whether this would be possible. – ali_m Feb 18 '15 at 13:10
• 1. See Edit. 2. vague and not a robust or particularly useful solution. 3. More information about your needs might help. Sorry, your question was a bit of a 2-parter and I think I was assuming that the only part you were having trouble with was the end part. removing the heartbeat. – Andrew Gallasch Feb 19 '15 at 13:20