7
$\begingroup$

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.

enter image description here

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:

enter image description here

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:

enter image description here

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: enter image description here

And the corresponding periodogram: enter image description here

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)

$\endgroup$
  • $\begingroup$ Starting from a series of movie frames, can you please elucidate more clearly how you are generating the PSD exactly? $\endgroup$ – Tarin Ziyaee Dec 16 '13 at 19:17
  • $\begingroup$ @user4619 very crudely - for each frame I've just computed the average pixel value to generate a vector x, then I take Fx = rfft(x), and get the power as abs(Fx)**2 $\endgroup$ – ali_m Dec 16 '13 at 19:29
  • $\begingroup$ You have a 2-D frame, and then you generate an average 1-D vector. Along x? Along y? $\endgroup$ – Tarin Ziyaee Dec 16 '13 at 19:35
  • $\begingroup$ @user4619 along both x and y - I reshape my movie into an nframes by npixels array, then average across all pixels $\endgroup$ – ali_m Dec 16 '13 at 19:38
  • $\begingroup$ Ok, thanks for that detail - it matters in the analysis. Please add this information to your post. $\endgroup$ – Tarin Ziyaee Dec 16 '13 at 19:54
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Isn't that essentially what I'm already doing with the notch filter? I still think it's not quite ideal, though, since this approach doesn't take into account information about the phase of the oscillation I'm trying to remove, which is constant over time. It seems to me like it should be possible to selectively remove the artifact without affecting 'genuine' signal that falls into the same frequency band. $\endgroup$ – ali_m Jun 24 '14 at 16:38
  • $\begingroup$ It's not exactly what you did: firstly, it's a much narrower notch; secondly, unlike the Butterworth filter, it causes no phase distortion. If you subtract the FFTs before and after modification, you get a single nonzero component with some amplitude, and the phase of the original spike. This is the signal we're removing in the time domain, i.e. a sinusoid of constant amplitude at exactly the frequency and phase of the spectral peak. If the underlying signal has energy in this region, it appears as noise in the estimate, but should wash out. $\endgroup$ – dpwe Jun 25 '14 at 0:19
  • $\begingroup$ I realize you may have meant enforcing the same phase across every pixel, so I added that to my answer. It's not magic, though - you can't separate the in-phase component of signal and noise, so you're always left with a residual 90\deg shifted from the estimated artifact. $\endgroup$ – dpwe Jun 25 '14 at 0:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.