Fourier-Raum Filterung

Question

Ich habe eine echte Vektor-Zeitreihe x der Länge T und ein Filter h der Länge t < < T. h ist ein Filter im Fourier-Raum, real und symmetrisch. Es ist ungefähr 1/f.

Ich möchte x mit h filtern, um y zu erhalten.

Angenommen t == T und FFTs der Länge T könnten in den Speicher passen (keine davon ist wahr). Zu meiner gefilterten x in Python, ich tun würde:

import numpy as np 
from scipy.signal import fft, ifft 

y = np.real(np.ifft(np.fft(x) * h))) 

Da die Bedingungen nicht halten, habe ich versucht, den folgenden Hack:

1. auswählen padding Größe P < t/2, wählen Sie eine Blockgröße B, so dass B + 2P eine gute FFT-Größe ist
2. Skalierung h über Spline-Interpolation zur Größe B + 2P> t (h_skaliert)
3. y = []; Loop:
   • Nehmen Block der Länge B + 2P von x (genannt x_b)
   • Perform y_b = IFFT (fft (x_b) * h_scaled)
   • Tropfen padding P von jeder Seite des y_b und Verketten mit y
   • Select nächste x_b mit zuletzt überlappend von P

Ist das eine gute Strategie? Wie wähle ich die Polsterung P auf eine gute Weise aus? Was ist der richtige Weg, dies zu tun? Ich weiß nicht viel Signalverarbeitung. Das ist eine gute Chance zu lernen.

Ich benutze cuFFT, um Dinge zu beschleunigen, so wäre es toll, wenn der Großteil der Operationen FFTs sind. Das eigentliche Problem ist 3D. Ich mache mir auch keine Gedanken über Artefakte aus einem akausalen Filter.

Danke, Paul.

Answer 1

Sie sind auf dem richtigen Weg. Die Technik heißt overlap-save processing. Ist t kurz genug, dass FFTs dieser Länge in den Speicher passen? Wenn ja, können Sie Ihre Blockgröße B so wählen, dass B > 2*min(length(x),length(h)) und für eine schnelle Transformation sorgt. Dann, wenn Sie verarbeiten, lassen Sie die erste Hälfte von y_b fallen, anstatt von beiden Enden zu fallen.

Um zu sehen, warum Sie die erste Hälfte fallen lassen, denken Sie daran, dass die spektrale Multiplikation die gleiche ist wie die zirkuläre Faltung im Zeitbereich. Das Konvertieren mit der Null-Padded h erzeugt seltsame Glitchy Transients in der ersten Hälfte des Ergebnisses, aber in der zweiten Hälfte sind alle Transienten verschwunden, da der Circular-Wrap-Punkt in x mit dem Null-Teil von h ausgerichtet ist. Es gibt eine gute Erklärung dafür, mit Bildern, in "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold.

Es ist wichtig, dass Ihr Zeitdomänenfilter mindestens an einem Ende auf 0 reduziert wird, so dass Sie keine Klingelartefakte erhalten. Sie erwähnen, dass h real in der Frequenzdomäne ist, was bedeutet, dass es alle 0 Phase hat. Gewöhnlich wird ein solches Signal nur in einer kreisförmigen Weise kontinuierlich sein, und wenn es als ein Filter verwendet wird, wird eine Verzerrung über das gesamte Frequenzband hinweg erzeugt. Eine einfache Möglichkeit, einen vernünftigen Filter zu erstellen, besteht darin, ihn in der Frequenzdomäne mit 0-Phase, inverser Transformation und Rotation zu entwerfen.Zum Beispiel:

def OneOverF(N): 
    import numpy as np 
    N2 = N/2; #N has to be even! 
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1))) 
    hf = 1/(2*np.pi*x/N2) 
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error 
    htrot = np.roll(ht, N2) 
    htwin = htrot * np.hamming(N) 
    return ht, htrot, htwin 

(Ich bin ziemlich neu in Python, lassen Sie es mich wissen, ob es ein besserer Weg, dies zu kodieren).

Wenn Sie die Frequenzgänge von ht vergleichen, htrot und htwin, sehen Sie die folgende (x-Achse ist die normalisierte Frequenz bis zu pi):

ht, an der Spitze, hat viele Welligkeit . Dies liegt an der Diskontinuität am Rand. htrot, in der Mitte, ist besser, aber immer noch Welligkeit. htwin ist schön und glatt, auf Kosten der Abflachung bei einer etwas höheren Frequenz. Beachten Sie, dass Sie die Länge des geradlinigen Abschnitts erweitern können, indem Sie einen größeren Wert für N verwenden.

Ich schrieb über das Diskontinuitätsproblem und schrieb auch ein Matlab/Octave-Beispiel in another SO question, wenn Sie weitere Details anzeigen möchten.