DSP - Filtern im Frequenzbereich über FFT

Question

Ich habe um ein wenig gespielt, mit der Exocortex Umsetzung der FFT, aber ich habe einige Probleme.

Wenn ich die Amplituden des Frequenzbins modifiziere, bevor die IFFT das resultierende Signal Aufruf enthält einige Klicken und Knacken, insbesondere bei niedrigeren Frequenzen in dem Signal vorhanden sind (wie Trommeln oder Bass). Dies ist jedoch nicht der Fall sein, wenn ich alle die Behälter mit dem gleichen Faktor dämpfen.

Lassen Sie mich stellen ein Beispiel des Ausgangspuffers eines 4-Probe FFT:

// Bin 0 (DC)
FFTOut[0] = 0.0000610351563
FFTOut[1] = 0.0

// Bin 1
FFTOut[2] = 0.000331878662
FFTOut[3] = 0.000629425049

// Bin 2
FFTOut[4] = -0.0000381469727
FFTOut[5] =  0.0

// Bin 3, this is the first and only negative frequency bin.
FFTOut[6] =  0.000331878662
FFTOut[7] = -0.000629425049

Der Ausgang besteht aus Paaren von Schwimmern, die jeweils die Real- und imaginay Teile eines einzigen Behälters darstellt. So ist 0 (Array-Indizes 0, 1) würde stellen die Real- und Imaginärteile der DC-Frequenz. Wie Sie sehen können, Behälter 1 und 3 beide haben die gleichen Werte, (außer dem Zeichen des Im-Teil), also denke ich bin 3 die erste negative Frequenz ist, und schließlich Indizes (4, 5) wäre der letzte positiv sein Frequenz ist.

Dann den Frequenz-Bins zu dämpfen 1 das ist, was ich tue:

// Attenuate the 'positive' bin
FFTOut[2] *= 0.5;
FFTOut[3] *= 0.5;

// Attenuate its corresponding negative bin.
FFTOut[6] *= 0.5;
FFTOut[7] *= 0.5;

Für die eigentlichen Tests, die ich bin mit einer 1024-Länge FFT und ich alle Proben immer zur Verfügung stellen, so dass keine 0-Polsterung benötigt wird.

// Attenuate
var halfSize = fftWindowLength / 2;
float leftFreq = 0f;
float rightFreq = 22050f; 
for( var c = 1; c < halfSize; c++ )
{
    var freq = c * (44100d / halfSize);

    // Calc. positive and negative frequency indexes.
    var k = c * 2;
    var nk = (fftWindowLength - c) * 2;

    // This kind of attenuation corresponds to a high-pass filter.
    // The attenuation at the transition band is linearly applied, could
    // this be the cause of the distortion of low frequencies?
    var attn = (freq < leftFreq) ? 
                    0 : 
                    (freq < rightFreq) ? 
                        ((freq - leftFreq) / (rightFreq - leftFreq)) :
                        1;

    // Attenuate positive and negative bins.
    mFFTOut[ k ] *= (float)attn;
    mFFTOut[ k + 1 ] *= (float)attn;
    mFFTOut[ nk ] *= (float)attn;
    mFFTOut[ nk + 1 ] *= (float)attn;
}

Offensichtlich etwas mache ich falsch, aber kann nicht herausfinden, was.

Ich mag nicht die FFT-Ausgabe als Mittel verwenden, um einen Satz von FIR-Koeffizienten zu erzeugen, da ich einen sehr dynamischen Equalizer zu implementieren bin versucht.

Was ist der richtige Weg, um Filter im Frequenzbereich? was ich bin fehlt?

Auch ist es wirklich so gut abzuschwächen negativen Frequenzen benötigt? Ich habe eine FFT-Implementierung gesehen, wo neg. Frequenzwerte werden vor der Synthese auf Null gesetzt.

Vielen Dank im Voraus.

Answer 1

Es gibt zwei Probleme: die Art und Weise Sie die FFT verwenden, und die besondere Filter

.

Die Filterung wird traditionell als Faltung im Zeitbereich implementiert. Du hast Recht, dass die Spektren der Eingangs- und Filtersignale multipliziert entspricht. Wenn Sie jedoch die diskrete Fourier-Transformation (DFT) verwenden (mit einer Fast-Fourier-Algorithmus für die Geschwindigkeit Transformation implementiert), die Sie tatsächlich eine abgetastete Version des wahren Spektrums berechnen. Dies hat viele Auswirkungen, aber die am meisten relevante Filterung ist die Implikation, dass das Zeitsignal periodisch ist.

Hier ist ein Beispiel. Betrachten wir ein sinusförmiges Eingangssignal mit x 1,5 Zyklen in der Zeit und ein einfacher Tiefpassfilter h. In Matlab / Octave Syntax:

N = 1024;
n = (1:N)'-1; %'# define the time index
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4
h = [h./sum(h)]; %# normalize DC gain

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra
y = real(y); %# due to numerical error, y has a tiny imaginary part
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here
plot(y);

Und hier ist die grafische Darstellung:

Der Glitch am Anfang des Blocks ist nicht das, was wir alle erwarten. Aber wenn Sie fft(x) betrachten, macht es Sinn. Die diskrete Fourier-Transformation nimmt das Signal innerhalb des Transformationsblockes periodisch ist. Soweit die DFT weiß, fragten wir nach einer Periode dies die Transformation:

Dies führt zu den ersten wichtigen Aspekt bei der Filterung mit DFTs: Sie implementieren tatsächlich Kreiskonvolution , nicht lineare Faltung. So ist das „Glitch“ in der ersten Graph ist nicht wirklich eine Störung, wenn Sie die Mathematik zu betrachten. So dann wird die Frage: Gibt es einen Weg, um die Periodizität zu umgehen? Die Antwort lautet ja: Verwendung Overlap-save Verarbeitung . Im Wesentlichen berechnen Sie N-Langprodukte wie oben, aber nur halten N / 2 Punkte.

Nproc = 512;
xproc = zeros(2*Nproc,1); %# initialize temp buffer
idx = 1:Nproc; %# initialize half-buffer index
ycorrect = zeros(2*Nproc,1); %# initialize destination
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer
    idx = idx + Nproc; %# step half-buffer index
end

Und hier ist die grafische Darstellung der ycorrect:

Dieses Bild macht Sinn - wir erwarten einen Anfahr-Übergang aus dem Filter, wird das Ergebnis setzt sich in den stationären Zustand Sinus Antwort. Beachten Sie, dass jetzt x beliebig lang sein können. Die Begrenzung ist Nproc > 2*min(length(x),length(h)).

Jetzt auf die zweite Frage: die jeweiligen Filter. In der Schleife, erstellen Sie einen Filter, der im Wesentlichen Spektrum ist ist H = [0 (1:511)/512 1 (511:-1:1)/512]'; Wenn Sie hraw = real(ifft(H)); plot(hraw) tun, erhalten Sie:

Es ist schwer zu sehen, aber es gibt eine Reihe von Nicht-Null-Punkte am äußersten linken Rand des Diagramms, und dann ein paar mehr am rechten Rand. Mit Octave eingebaute in freqz Funktion bei dem Frequenzgang suchen wir (indem freqz(hraw)) finden Sie unter:

Die Größe Reaktion hat viele Wellen aus dem Hochpass Hüllkurve auf Null. Auch hier ist die Periodizität, die mit der DFT bei der Arbeit. Soweit die DFT betrifft, hraw wiederholt immer und immer wieder. Aber wenn Sie nehmen eine Periode von hraw, wie freqz tut, sein Spektrum ist ganz anders als die periodische Version der.

Lassen Sie sich also ein neues Signal zu definieren: hrot = [hraw(513:end) ; hraw(1:512)]; Wir drehen einfach die rohe DFT Ausgabe, die es kontinuierlich innerhalb des Blockes zu machen. Nun ist aussehen läßt bei dem Frequenzgang mit freqz(hrot):

Viel besser. Die gewünschte Hüllkurve ist es, ohne alle Wellen. Natürlich ist die Umsetzung jetzt nicht so einfach ist, haben Sie eine vollständige komplexe Multiplikation mit fft(hrot) zu tun, anstatt nur jede komplexe bin Skalierung, aber zumindest werden Sie die richtige Antwort.

Beachten Sie, dass für die Geschwindigkeit, dann würden Sie in der Regel die DFT des gepolsterten h vorab berechnen, verließ ich es allein in der Schleife mehr easily vergleicht mit dem Original.

Answer 2

Ihr Hauptproblem ist, dass die Frequenzen sind nicht gut kurze Zeitintervalle definiert über . Dies gilt insbesondere für niedrigen Frequenzen, weshalb Sie das Problem der meisten es bemerken.

Wenn Sie also wirklich nehmen kurze Segmente aus dem Klang Zug, und dann diese Sie filtern, die gefilterte Segmente wont Filter in einer Weise, die eine kontinuierliche Wellenform erzeugt, und man hört die zwischen den Segmenten springt und das ist, was die erzeugt Klicks, die Sie hier.

Zum Beispiel, einige vernünftige Zahlen nehmen: ich mit einer Wellenform bei 27,5 Hz (A0 auf einem Klavier) beginnen bei 44.100 Hz digitalisiert, wird es wie folgt aussehen (wo der rote Teil 1024 Abtastungen lang ist):

Also zuerst werden wir mit einem Tiefpass von 40Hz starten. Also, da die ursprüngliche Frequenz kleiner als 40 Hz, ein Tiefpassfilter mit einem 40Hz cut-off soll eigentlich keine Wirkung hat, und wir werden eine Ausgabe erhalten, dass fast genau die Eingabe übereinstimmt. Richtig? Falsch, falsch, falsch - und das ist im Grunde der Kern des Problems. Das Problem ist, dass für die kurzen Abschnitte der Idee von 27,5 Hz ist nicht klar definiert, und kann nicht gut in der DFT dargestellt werden.

Das 27,5 Hz ist nicht besonders sinnvoll in dem kurzen Segment unten, indem man die DFT in der Figur zu sehen ist. Beachten Sie, dass obwohl die längere Segment DFT (schwarze Punkte) zeigt einen Peak bei 27,5 Hz, die kurze (rote Punkte) nicht.

ist klar, dann unter 40Hz Filterung, erfassen wird nur die DC-Offset, und das Ergebnis des 40Hz Tiefpassfilter ist in grün unten.

Die blaue Kurve (aufgenommen mit einem 200 Hz cut-off) beginnt, viel besser zusammenpassen. Aber beachten Sie, dass es nicht die niedrigen Frequenzen, die sie zusammenpassen gut, aber die Aufnahme von hohen Frequenzen machen. Es ist nicht, bis wir jede Frequenz möglich in dem kurzen Segment bis zu 22 kHz sind, dass wir endlich eine gute Darstellung der ursprünglichen Sinuswelle erhalten.

Der Grund für all dies ist, dass ein kleines Segment einer 27,5 Hz Sinuswelle ist nicht eine 27,5 Hz Sinuswelle, und es ist DFT nicht viel mit 27,5 Hz zu tun.

Answer 3

mildernden Sie den Wert der DC Frequenz Probe auf Null? Es scheint, dass Sie nicht in Ihrem Beispiel überhaupt sind zu dämpfen. Da Sie einen Hochpassfilter implementieren, müssen Sie auch den DC-Wert auf Null setzen.

Dies würde Niederfrequenzverzerrung erklären. Sie würden bei niedrigeren Frequenzen im Frequenzgang eine Menge Welligkeit, wenn die DC-Wert nicht Null ist wegen des großen Übergangs ist.

Hier ist ein Beispiel in MATLAB / Octave zu zeigen, was passiert sein könnte:

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,4) linspace(1,0,8) zeros(1,3) 1 zeros(1,4) linspace(0,1,8) ones(1,4)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude');

Beachten Sie, dass in meinem Code, ich ein Beispiel für den DC-Wert erschaffe nicht Null sein, gefolgt von einer abrupten Änderung auf Null und dann einer Rampe nach oben. Ich nehme dann die IFFT in den Zeitbereich zu transformieren. Dann fahre ich eine Null-gepolsterte fft (die automatisch von MATLAB erfolgt, wenn Sie in einer FFT-Größe größer als das Eingangssignal übergeben) auf dem Zeitbereichssignal. Das Zero-Padding in der Zeitdomäne führt Interpolation in der Frequenzdomäne. Mit diesem können wir sehen, wie die Filter zwischen Filterproben reagieren werden.

Eines der wichtigsten Dinge zu erinnern ist, dass, obwohl Sie Filterantwortwerte bei bestimmten Frequenzen setzen, indem sie die Ausgänge der DFT Abschwächen, diese Garantien nichts für Frequenzen zwischen Abtastpunkten auftreten. Dies bedeutet, die abrupte Änderung, desto mehr Überschwingen und Pendel zwischen Proben auftreten.

Jetzt Ihre Frage zu beantworten, wie diese Filterung getan werden soll. Es gibt eine Reihe von Möglichkeiten, aber eine der am einfachsten zu implementieren und zu verstehen, wird das Fenster Design-Methode. Das Problem mit dem aktuellen Design ist, dass die Übergangsbreite ist riesig. Die meiste Zeit, werden Sie so schnell wie möglich von Übergängen, mit möglichst wenig Welligkeit wie möglich.

Im nächsten Code, ich werde einen idealen Filter erstellen und die Antwort angezeigt werden:

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,8) zeros(1,16) ones(1,8)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude'); 

Beachten Sie, dass es eine Menge Schwingung durch die abrupten Veränderungen verursacht.

Die FFT oder diskrete Fourier-Transformation ist eine abgetastete Version der Fourier-Transformation. Die Fourier-Transformation wird auf ein Signal über den fortlaufenden Bereich -Infinity bis unendlich aufgebracht, während die DFT über eine endliche Anzahl von Proben angewandt werden. Dieser Effekt führt zu einem Quadrat Windowing (Trunkierung) in der Zeitdomäne, wenn die DFT verwenden, da wir nur mit einer begrenzten Anzahl von Proben zu tun sind. Leider ist die DFT einer Rechteckwelle eine sinc Typ Funktion (sin (x) / x).

Das Problem bei mit scharfen Übergängen im Filter (schneller Sprung von 0 auf 1 in einer Probe) ist, dass dies eine sehr lange Antwort im Zeitbereich hat, der durch ein quadratisches Fenster abgeschnitten wird. So, um Hilfe zu minimieren dieses Problem können wir das Zeitsignal durch ein allmähliches Fenster multiplizieren. Wenn wir mehrfach ein Hanning-Fenster, indem Sie die Zeile:

x = x .* hanning(1,N).';

nach der IFFT nehmen wir diese Antwort bekommen:

Also ich würde empfehlen, versucht das Fenster Design-Methode zu implementieren, da es relativ einfach ist (es gibt bessere Möglichkeiten, aber sie noch komplizierter). Da Sie einen Equalizer implementieren, gehe ich davon aus Sie die Dämpfungen on the fly zu ändern, um die Lage sein wollen, so würde ich vorschlagen, die Berechnung und die Filter in der Frequenzdomäne zu speichern, wenn es in den Parametern eine Änderung ist, und dann kann man es einfach anwenden an jedes Eingangsaudio durch nehmen der FFT des Eingangspuffers, Multiplizieren mit Ihrem Frequenzbereichsfilterproben-Puffer und dann wird die DurchführungIFFT in den Zeitbereich zurück zu bekommen. Dies wird viel effizienter als alle für jede Probe der Verzweigung ist Sie tun.

Answer 4

Zuerst über die Normalisierung: das ein bekanntes (nicht) Problem. Die DFT / IDFT würde einen Faktor erfordern 1 / sqrt (N) (abgesehen von der Standard-Cosinus / Sinus-Faktoren) in jedem (eine inverse direkt), um sie simmetric und wirklich umkehrbar zu machen. Eine weitere Möglichkeit ist eine von ihnen zu teilen (die direkte oder inverse) von N , das ist eine Sache der Bequemlichkeit und Geschmack. Oft sind die FFT-Routinen diese Normalisierung nicht durchführen, wird der Benutzer erwartet von ihm und normalisieren sich bewusst sein, wie er es vorzieht. Siehe

Zweitens: in einer (sagen wir) 16-Punkt-DFT, was Sie rufen die bin 0 entspräche der Nullfrequenz (DC), bin 1 niedrige Freq .. . ist 4 Medium freq, ist 8 zu der höchsten Frequenz und bins 9 ... 15 zu den "negativen Frequenzen". In Ihrem Beispiel dann, ist 1 ist eigentlich sowohl die niedrige Frequenz und mittlere Frequenz. Abgesehen von dieser Überlegung ist es nichts konzeptionell falsch in Ihrer „Entzerrung“. Ich verstehe nicht, was Sie unter „das Signal bei niedrigeren Frequenzen verzerrt wird“ . Wie beobachten Sie das?