Entfernen eines sinusförmigen Artefakts aus einer Reihe von Filmbildern

Ich mache eine Post-hoc-Analyse eines Datensatzes, der aus einer Reihe von Filmbildern besteht, die durch ein stark periodisches Artefakt kontaminiert sind. Ich möchte dieses Artefakt aus meinen Frames entfernen.

Um das Plotten zu vereinfachen, habe ich gerade mein Array Mvon Pixelwerten auf umgeformt [nframes, npixels]und dann über alle Pixelwerte gemittelt, um einen 1D-Vektor zu erhalten m. So sieht dieses Signal im Zeitbereich aus. Sie können die Schwingung im gezoomten Einschub ziemlich deutlich sehen.

Ich machte dann ein Periodogramm durch Aufnehmen Fm = rfft(m)und zeichnete abs(Fm)**2gegen die Frequenz. Ich sehe eine sehr scharfe Spitze bei ~ 1,5 Hz:

Neben der zeitlichen Periodizität scheint dieses Artefakt auch eine schwächere räumliche Komponente zu haben, da beim exakten Spitzenfrequenzwert eine gleichmäßige Phasenänderung über die x-Achse meiner Frames zu erfolgen scheint, so dass Pixel auf der rechts neigen dazu, Pixel links hinterherzuhinken:

Als Brute-Force-Ansatz habe ich versucht, jedes Pixel im Zeitbereich mit einem auf 1,5 Hz zentrierten Sperrfilter zu filtern. Ich habe einen Butterworth-Filter der Ordnung 4 mit kritischen Frequenzen von 1,46 und 1,52 Hz verwendet (ich bin mit dem Filterdesign nicht vertraut, daher bin ich mir sicher, dass es geeignetere Optionen gibt).

So sieht das mittlere Pixelsignal nach dem Filtern aus:

Und das entsprechende Periodogramm:

Der Sperrfilter reduziert das Artefakt einigermaßen gut, aber da es im Grunde genommen wie eine reine stationäre Sinuskurve aussieht, kann ich nicht anders, als zu glauben, dass ich es besser machen könnte, als nur diesen Teil des Frequenzraums zu dämpfen.

Meine anfängliche (sehr naive) Idee war, etwas zu tun wie:

1. Holen Sie sich die Frequenz, Phase und Amplitude der Schwingung aus dem Fourier-Spektrum für jedes Pixel im Film
2. Rekonstruieren Sie die Schwingung im Zeitbereich
3. Subtrahieren Sie es von den Filmbildern

Mir ist klar, dass dies normalerweise nicht der Fall ist, da Interferenzen normalerweise nicht so spektral rein und zeitlich stationär sind, aber ich frage mich, ob dies in meinem Fall sinnvoll sein könnte.

Daten

Voller 16-Bit-TIFF-Stapel (~ 2 GB unkomprimiert)

Räumlich dezimierte 8-Bit-Version (~ 35 MB unkomprimiert)

Ihre vorgeschlagene Lösung - Berechnen einer Sinuskurve im Zeitbereich basierend auf dem Peak in der FFT und anschließendes Subtrahieren - sollte funktionieren, aber es gibt einen einfacheren Weg, im Wesentlichen dasselbe zu tun: Ändern Sie diesen Peakwert in der FFT und nehmen Sie dann die Umkehrung verwandeln.

Für Ihr gerastertes Video M[nframes, npixels]finden Sie also das Frequenzfach, in dem sich das Artefakt befindet, und reduzieren es dann systematisch (z. B. stellen Sie seine Größe auf den Durchschnitt seiner Nachbarn ein) für jedes 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)

Dies sollte funktionieren, wenn das Artefakt eine genau konstante Amplitude und Frequenz aufweist und seine Frequenz direkt auf ein Submultiplikator der Sequenzlänge (dh die durch die FFT dargestellten Sinuskurven) fällt. Im Allgemeinen möchten Sie möglicherweise ein oder zwei Bins auf beiden Seiten abflachen, badbinum einen etwas breiteren Satz von Schmalbandverfälschungen zu behandeln, z

# ...

   # 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])
   # ...

Wenn Sie die von jedem Pixel entfernte Komponente so einschränken möchten, dass dieselbe Frequenz und Phase des Artefakts in der mittleren Intensität erfasst wird, können Sie nur die Projektion der badbinGröße auf diese Phase entfernen , z

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)
   # ...

Es ist zu beachten, dass die resultierende Komponente bei jedem Pixel badbinjetzt immer um 90 ° phasenverschoben (orthogonal) von der globalen badbinphaseKomponente ist - jede Signalkomponente mit genau dieser Frequenz und Phase kann nicht vom Artefakt getrennt werden.