Inverser Kurzzeit-Fourier-Transformations-Algorithmus in Worten beschrieben

Ich versuche begrifflich zu verstehen, was passiert, wenn die vorwärts und rückwärts gerichteten Kurzzeit-Fourier-Transformationen (STFT) auf ein diskretes Zeitdomänensignal angewendet werden. Ich habe das klassische Papier von Allen und Rabiner ( 1977 ) sowie einen Wikipedia-Artikel ( Link ) gefunden. Ich glaube, dass es hier noch einen weiteren guten Artikel gibt .

Ich interessiere mich für die Berechnung der Gabor-Transformation, die nichts anderes als die STFT mit einem Gaußschen Fenster ist.

Folgendes verstehe ich über die vorwärts gerichtete STFT:

1. Aus dem Signal wird eine Teilsequenz ausgewählt, die aus Zeitbereichselementen besteht.
2. Die Teilsequenz wird mit einer Fensterfunktion unter Verwendung einer Punkt-für-Punkt-Multiplikation im Zeitbereich multipliziert.
3. Die multiplizierte Teilsequenz wird mit der FFT in den Frequenzbereich übernommen.
4. Durch Auswahl aufeinanderfolgender überlappender Teilsequenzen und Wiederholung der obigen Prozedur erhalten wir eine Matrix mit m Zeilen und n Spalten. Jede Spalte ist die zu einem bestimmten Zeitpunkt berechnete Teilsequenz. Dies kann verwendet werden, um ein Spektrogramm zu berechnen.

Für die inverse STFT sprechen die Arbeiten jedoch von einer Summierung über die überlappenden Analyseabschnitte. Ich finde es sehr schwierig zu visualisieren, was hier wirklich vor sich geht. Was muss ich tun, um die inverse STFT berechnen zu können (in schrittweiser Reihenfolge wie oben)?

STFT weiterleiten

Ich habe eine Zeichnung erstellt, die zeigt, was meiner Meinung nach für die vorwärtsgerichtete STFT vor sich geht. Was ich nicht verstehe, ist, wie man die einzelnen Teilsequenzen zusammensetzt, damit ich die ursprüngliche Zeitsequenz zurückerhalte. Könnte jemand diese Zeichnung modifizieren oder eine Gleichung angeben, die zeigt, wie die Teilsequenzen hinzugefügt werden?

Inverse Transformation

Folgendes verstehe ich über die inverse Transformation. Jedes nachfolgende Fenster wird mit der IFFT in den Zeitbereich zurückgeführt. Dann wird jedes Fenster um die Schrittgröße verschoben und zum Ergebnis der vorherigen Verschiebung hinzugefügt. Das folgende Diagramm zeigt diesen Vorgang. Die summierte Ausgabe ist das Zeitdomänensignal.

Code-Beispiel

Der folgende Matlab-Code generiert ein synthetisches Zeitdomänensignal und testet dann den STFT-Prozess, um zu demonstrieren, dass das Inverse das Dual der Vorwärtstransformation innerhalb eines numerischen Rundungsfehlers ist. Anfang und Ende des Signals sind mit Nullen aufgefüllt, um sicherzustellen, dass die Mitte des Fensters am ersten und letzten Element des Zeitbereichssignals liegt.

$N + N_0 - 1$$N_0$

% The code computes the STFT (Gabor transform) with step size = 1
% This is most useful when modifications of the signal is required in
% the frequency domain

% The Gabor transform is a STFT with a Gaussian window (w_t in the code)

% written by Nicholas Kinar

% Reference:
% [1] J. B. Allen and L. R. Rabiner, 
% “A unified approach to short-time Fourier analysis and synthesis,” 
% Proceedings of the IEEE, vol. 65, no. 11, pp. 1558 – 1564, Nov. 1977.

% generate the signal
mm = 8192;                  % signal points
t = linspace(0,1,mm);       % time axis

dt = t(2) - t(1);           % timestep t
wSize = 101;                % window size


% generate time-domain test function
% See pg. 156
% J. S. Walker, A Primer on Wavelets and Their Scientific Applications, 
% 2nd ed., Updated and fully rev. Boca Raton: Chapman & Hall/CRC, 2008.
% http://www.uwec.edu/walkerjs/primer/Ch5extract.pdf
term1 = exp(-400 .* (t - 0.2).^2);
term2 = sin(1024 .* pi .* t);
term3 = exp(-400.*(t- 0.5).^2);
term4 = cos(2048 .* pi .* t);
term5 = exp(-400 .* (t-0.7).^2);
term6 = sin(512.*pi.*t) - cos(3072.*pi.*t);
u = term1.*term2  + term3.*term4 + term5.*term6; % time domain signal
u = u';

figure;
plot(u)

Nmid = (wSize - 1) / 2 + 1;    % midway point in the window
hN = Nmid - 1;                 % number on each side of center point       


% stores the output of the Gabor transform in the frequency domain
% each column is the FFT output
Umat = zeros(wSize, mm);     


% generate the Gaussian window 
% [1] Y. Wang, Seismic inverse Q filtering. Blackwell Pub., 2008.
% pg. 123.
T = dt * hN;                    % half-width
sp = linspace(dt, T, hN); 
targ = [-sp(end:-1:1) 0 sp];    % this is t - tau
term1 = -((2 .* targ) ./ T).^2;
term2 = exp(term1);
term3 = 2 / (T * sqrt(pi));
w_t = term3 .* term2;
wt_sum = sum ( w_t ); % sum of the wavelet


% sliding window code
% NOTE that the beginning and end of the sequence
% are padded with zeros 
for Ntau = 1:mm

    % case #1: pad the beginning with zeros
    if( Ntau <= Nmid )
        diff = Nmid - Ntau;
        u_sub = [zeros(diff,1); u(1:hN+Ntau)];
    end

    % case #2: simply extract the window in the middle
    if (Ntau < mm-hN+1 && Ntau > Nmid)
        u_sub = u(Ntau-hN:Ntau+hN);
    end

    % case #3: less than the end
    if(Ntau >= mm-hN+1)
        diff = mm - Ntau;
        adiff = hN - diff;
        u_sub = [ u(Ntau-hN:Ntau+diff);  zeros(adiff,1)]; 
    end   

    % windowed trace segment
    % multiplication in time domain with
    % Gaussian window  function
    u_tau_omega = u_sub .* w_t';

    % segment in Fourier domain
    % NOTE that this must be padded to prevent
    % circular convolution if some sort of multiplication
    % occurs in the frequency domain
    U = fft( u_tau_omega );

    % make an assignment to each trace
    % in the output matrix
    Umat(:,Ntau) = U;

end

% By here, Umat contains the STFT (Gabor transform)

% Notice how the Fourier transform is symmetrical 
% (we only need the first N/2+1
% points, but I've plotted the full transform here
figure;
imagesc( (abs(Umat)).^2 )


% now let's try to get back the original signal from the transformed
% signal

% use IFFT on matrix along the cols
us = zeros(wSize,mm);
for i = 1:mm 
    us(:,i) = ifft(Umat(:,i));
end

figure;
imagesc( us );

% create a vector that is the same size as the original signal,
% but allows for the zero padding at the beginning and the end of the time
% domain sequence
Nuu = hN + mm + hN;
uu = zeros(1, Nuu);

% add each one of the windows to each other, progressively shifting the
% sequence forward 
cc = 1; 
for i = 1:mm
   uu(cc:cc+wSize-1) = us(:,i) + uu(cc:cc+wSize-1)';
   cc = cc + 1;
end

% trim the beginning and end of uu 
% NOTE that this could probably be done in a more efficient manner
% but it is easiest to do here

% Divide by the sum of the window 
% see Equation 4.4 of paper by Allen and Rabiner (1977)
% We don't need to divide by L, the FFT transform size since 
% Matlab has already taken care of it 
uu2 = uu(hN+1:end-hN) ./ (wt_sum); 

figure;
plot(uu2)

% Compare the differences bewteen the original and the reconstructed
% signals.  There will be some small difference due to round-off error
% since floating point numbers are not exact
dd = u - uu2';

figure;
plot(dd);

Das STFT-Transformationspaar kann durch 4 verschiedene Parameter charakterisiert werden:

1. FFT-Größe (N)
2. Schrittweite (M)
3. Analysefenster (Größe N)
4. Synthesefenster (Größe N)

Der Prozess ist wie folgt:

1. Nehmen Sie N (fft size) Samples vom aktuellen Eingangsort auf
2. Analysefenster anwenden
3. Mach die FFT
4. Tun Sie alles, was Sie im Frequenzbereich tun möchten
5. Inverse FFT
6. Synthesefenster anwenden
7. Zur Ausgabe am aktuellen Ausgabeort hinzufügen
8. Erweitern Sie die Eingabe- und Ausgabeposition um M (Schrittgröße) Samples

Der Überlappungsadditionsalgorithmus ist dafür ein gutes Beispiel. In diesem Fall ist die Schrittgröße N, die FFT-Größe 2 · N, das Analysefenster ist rechteckig mit N Einsen, gefolgt von N Nullen, und die Synthesefenster sind einfach alle Einsen.

Dafür gibt es viele andere Möglichkeiten und unter bestimmten Bedingungen wird die Vorwärts- / Rückwärtsübertragung vollständig rekonstruiert (dh Sie können das ursprüngliche Signal zurückerhalten).

Der entscheidende Punkt hierbei ist, dass jede Ausgangsstichprobe typischerweise additive Beiträge von mehr als einer inversen FFT erhält. Die Ausgabe muss über mehrere Frames hinweg akkumuliert werden. Die Anzahl der beitragenden Frames ergibt sich einfach aus der FFT-Größe geteilt durch die Schrittgröße (ggf. aufgerundet).

Sieben Jahre nachdem diese Frage zum ersten Mal gestellt wurde, bin ich in diese Verwirrung geraten, ähnlich wie bei Nicholas Kinar. Hier möchte ich einige "inoffizielle" und "nicht vollständig gesicherte" persönliche Wahrnehmungsideen und Erklärungen geben.

Der Titel der folgenden Aussagen ist zur besseren Verständlichkeit übertrieben.

1. Der Vorwärtsprozess von STFT ist nicht wirklich dazu gedacht, das ursprüngliche Signal zu erhalten.
   • Bei Verwendung von STFT mit einem nicht trivialen Fenster (nicht alle) ist das Eingangssignal für FFT eine verzerrte / gedehnte Version des ursprünglichen Signalfragments.
   • Dies ist gut für die Feature-Extraktion, bei der nutzlose / redundante Daten herausgefiltert werden. Wie bei der Silbenerkennung sind nicht alle zeitlichen Daten erforderlich, um bestimmte Töne in einer Sprache zu erkennen.
   • Die Spitze im Fenstervektor stellt die Minderheit der Positionen in einem Audiosignal dar, auf die Algorithmen achten sollten.
2. Das rohe Ergebnis der inversen STFT könnte also etwas sein, das wir möglicherweise nicht intuitiv erwarten.
   • Es sollten die mit Fenstern versehenen Signalfragmente sein, nach denen der Ifft von STFT-Funktionen aussehen sollte.
3. Um die ursprünglichen Signalfragmente ohne Fenster zu erhalten, kann ein Inversfenster auf die Rohausgabe von ifft angewendet werden.
   • Es ist einfach, eine Zuordnungsfunktion zu entwerfen, mit der der Effekt des Hann / Hamming-Fensters rückgängig gemacht werden kann.
4. Das Synthesefenster ist dann beteiligt, um zeitliche Fragmentierungsüberlappungen zu behandeln
   • Da die ursprünglichen Signalfragmente ohne Fenster als bereits erhalten angesehen werden können, können beliebige "Übergangsgewichtungen" verwendet werden, um die überlappenden Teile zu interpolieren.
5. Wenn Sie bedenken möchten, dass das FFT einer Fenstersprache die schwachen Signale möglicherweise weniger beachtet, aber diese starken Signale verehrt, gibt es möglicherweise eine Möglichkeit, entsprechende Synthesefenster zu entwerfen.
6. Ein einfacher Algorithmus zur Erzeugung eines Synthesefensters kann auch unter Anwendung der folgenden Prinzipien angegeben werden:
   • Gewichten Sie die Positionen im Synthesefenster höher, wenn der Wert des Analysefensters für diese Position hoch ist, verglichen mit anderen Fragmenten, die diese Position überlappen.
   • Verringern Sie die Positionen im Synthesefenster, wenn der Wert des Analysefensters für diese Position niedrig ist, und andere überlappende Fragmente berücksichtigen diese Position mit einem größeren Wert des Analysefensters stärker.