Removing a sinusoidal artifact from a set of movie frames

Question

I am doing some post-hoc analysis of a dataset consisting of a series of movie frames that are contaminated by a strongly periodic artifact. I would like to remove this artifact from my frames.

For ease of plotting I've just reshaped my array M of pixel values to [nframes, npixels], then averaged over all pixel values to give me a 1D vector m. Here's what this signal looks like in the time domain. You can see the oscillation quite clearly in the zoomed inset.

I then made a periodogram by taking Fm = rfft(m), and plotted abs(Fm)**2 against frequency. I see a very sharp peak at ~1.5Hz:

As well as the temporal periodicity, there also seems to be a weaker spatial component to this artifact, since at the exact peak frequency value there seems to be a smooth variation in phase across the x-axis of my frames, so that pixels on the right tend to lag pixels on the left:

As a brute force approach, I've tried just filtering each pixel in the time domain using a notch filter centred on 1.5Hz. I used an order 4 Butterworth filter with critical frequencies 1.46 and 1.52Hz (I'm not well versed in filter design, so I'm sure there may be more appropriate choices).

Here's what the mean pixel signal looks like after filtering:

And the corresponding periodogram:

The notch filter does a reasonably good job of reducing the artifact, but since it basically looks like a pure stationary sinusoid I can't help but think that I could do better than just attenuating that part of frequency space.

My initial (very naive) idea was to do something like:

1. Get the frequency, phase and amplitude of the oscillation from the Fourier spectrum for each pixel in the movie
2. Reconstruct the oscillation in the time domain
3. Subtract it from the movie frames

I realise this isn't something people usually do, since interference usually isn't so spectrally pure and temporally stationary, but I wonder if it might make sense in my case?

Data

Full 16bit TIFF stack (~2GB uncompressed)

Spatially decimated 8bit version (~35MB uncompressed)

Answer 1

Your proposed solution - calculating a sinusoid in the time domain based on the peak in the FFT, then subtracting it - should work, but there's an easier way to do essentially the same thing: modify that peak value in the FFT, then take the inverse transform.

So, for your rasterized video M[nframes, npixels] you find the frequency bin holding the artifact, then systematically flatten it (e.g., set its magnitude to the average of its neighbors) for every pixel:

import numpy as np
nframes, npixels = np.shape(M)
# Identify the bin containing the sinusoidal artifact
# Use the average intensity for each image
m = np.mean(M, axis=1)
# Calculate the FFT
Fm = np.fft.rfft(m)
# Find the largest bin away from the low-frequency region
lowfreq = 100  # or something
badbin = lowfreq + np.argmax(Fm[lowfreq:]**2)

# Now adjust the amplitude of that bin in the FFT of each pixel
for pixel in range(npixels):
   Fpix = np.fft.rfft(M[:, pixel])
   # Scale magnitude of artifact bin to be the mean of its neighbors
   Fpix[badbin] *= np.mean(np.absolute(Fpix[[badbin-1, badbin+1]]))/np.absolute(Fpix[badbin])
   # Rewrite the time sequence of that pixel
   M[:,pixel] = np.fft.irfft(Fpix)

This should work if the artifact is exactly constant amplitude and frequency, and its frequency falls right on to a submultiple of the sequence length (i.e., the sinusoids represented by the FFT). In general, you might want to flatten out one or two bins either side of badbin to deal with a slightly broader set of narrowband corruptions, e.g.

# ...

   # Scale magnitude of artifact binS to be the mean of neighbors
   spread = 3  # flatten bins from (badbin - (spread-1)) to (badbin + (spread-1))
   # target value for new bins
   targetmag = np.mean(np.absolute(Fpix[[badbin-spread, badbin+spread]]))
   bins = range(badbin - (spread-1), badbin + spread)
   Fpix[bins] *= targetmag/np.abs(Fpix[bins])
   # ...

If you want to constrain the component removed from each pixel to have the same frequency and phase of the artifact detected in the mean intensity, you could remove just the projection of the badbin magnitude onto that phase, e.g.

badbinphase = np.angle(Fm[badbin])
# ...

   Ncomponent = np.abs(Fpix[badbin])*np.cos(np.angle(Fpix[badbin]) - badbinphase)
   Fpix[badbin] -= Ncomponent * np.exp(0+1j * badbinphase)
   # ...

Note that the resulting component at badbin will now always be 90\deg phase shifted (orthogonal) from the global badbinphase in every pixel - any signal component at exactly that frequency and phase cannot be separated from the artifact.