Algorytm odwrotnej krótkotrwałej transformacji Fouriera opisany słowami

Próbuję pojęciowo zrozumieć, co się dzieje, gdy do dyskretnego sygnału w dziedzinie czasu stosowane są odwrotne transformaty Fouriera (STFT) do przodu i do tyłu. Znalazłem klasyczny artykuł Allena i Rabinera ( 1977 ), a także artykuł w Wikipedii ( link ). Uważam, że można tu znaleźć jeszcze jeden dobry artykuł .

Interesuje mnie obliczenie transformacji Gabora, która jest niczym innym jak STFT z oknem Gaussa.

Oto, co rozumiem na temat przekazywania STFT:

1. Podsekwencja jest wybierana z sygnału, złożonego z elementów w dziedzinie czasu.
2. Podsekwencja jest mnożona przez funkcję okna za pomocą mnożenia punkt po punkcie w dziedzinie czasu.
3. Mnożona podsekwencja jest pobierana do dziedziny częstotliwości za pomocą FFT.
4. Wybierając kolejne nakładające się podsekwencje i powtarzając powyższą procedurę, otrzymujemy macierz z m wierszami i n kolumnami. Każda kolumna jest podsekwencją obliczoną w danym momencie. Można to wykorzystać do obliczenia spektrogramu.

Jednak w przypadku odwrotnego STFT, artykuły mówią o podsumowaniu na nakładających się sekcjach analizy. Bardzo trudno jest mi wyobrazić sobie, co naprawdę się tutaj dzieje. Co muszę zrobić, aby móc obliczyć odwrotny STFT (krok po kroku, jak wyżej)?

Naprzód STFT

Stworzyłem rysunek pokazujący, co według mnie dzieje się w przypadku STFT. Nie rozumiem, jak złożyć każdą z podsekwencji, aby odzyskać pierwotną sekwencję czasową. Czy ktoś może zmodyfikować ten rysunek lub podać równanie pokazujące, w jaki sposób dodawane są podsekwencje?

Odwrotna transformacja

Oto, co rozumiem na temat transformacji odwrotnej. Każde kolejne okno jest przenoszone z powrotem do dziedziny czasu za pomocą IFFT. Następnie każde okno jest przesuwane o rozmiar kroku i dodawane do wyniku poprzedniego przesunięcia. Poniższy schemat pokazuje ten proces. Sumowany sygnał wyjściowy jest sygnałem w dziedzinie czasu.

Przykład kodu

Poniższy kod Matlab generuje syntetyczny sygnał w dziedzinie czasu, a następnie testuje proces STFT, wykazując, że odwrotność jest podwójną transformacją do przodu, w ramach błędu numerycznego zaokrąglenia. Początek i koniec sygnału są zerowane, aby zapewnić, że środek okna może znajdować się na pierwszym i ostatnim elemencie sygnału w dziedzinie czasu.

UWAGA: zgodnie z Allen i Rabiner (1977), jeśli w dziedzinie częstotliwości zachodzi mnożenie w celu zmiany odpowiedzi częstotliwościowej, długość okna analizy musi być równa lub większa niż punktów, gdzie jest odpowiedzią filtra . Długość jest przedłużana przez wypełnienie zerowe. Kod testowy pokazuje po prostu, że odwrotność jest podwójną transformacją do przodu. Długość musi zostać przedłużona, aby uniknąć kołowego splotu.N $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);

Para transformacji STFT może charakteryzować się 4 różnymi parametrami:

1. Rozmiar FFT (N)
2. Rozmiar kroku (M)
3. Okno analizy (rozmiar N)
4. Okno syntezy (rozmiar N)

Proces przebiega następująco:

1. Pobierz próbki N (rozmiar fft) z bieżącej lokalizacji wejściowej
2. Zastosuj okno analizy
3. Wykonaj FFT
4. Rób wszystko, co chcesz robić w dziedzinie częstotliwości
5. Odwrotna FFT
6. Zastosuj okno syntezy
7. Dodaj do wyjścia w bieżącym miejscu wyjścia
8. Zwiększ lokalizację wejścia i wyjścia według próbek M (wielkość kroku)

Algorytm dodawania nakładki jest tego dobrym przykładem. W tym przypadku rozmiar kroku to N, rozmiar FFT to 2 * N, okno analizy jest prostokątne z N jedynymi, po których następuje N zer, a okna syntezy to po prostu wszystkie.

Jest na to wiele innych możliwości i pod pewnymi warunkami transfer do przodu / do tyłu jest w pełni odtwarzany (tzn. Można odzyskać oryginalny sygnał).

Kluczową kwestią jest to, że każda próbka wyjściowa zazwyczaj otrzymuje addytywny wkład z więcej niż jednego odwrotnego FFT. Dane wyjściowe muszą być gromadzone w wielu ramkach. Liczba ramek przyczyniających się jest po prostu podana przez rozmiar FFT podzielony przez rozmiar kroku (w razie potrzeby zaokrąglony w górę).

Siedem lat po tym, jak pytanie zostało postawione po raz pierwszy, wpadłem w zamieszanie podobne do @Nicholas Kinar. W tym miejscu chciałbym przedstawić kilka „nieoficjalnych” i „nie do końca zapewnionych poprawności” osobistych percepcyjnych pomysłów i wyjaśnień.

Tytuły poniższych stwierdzeń są przesadzone w celu zwiększenia zrozumiałości.

1. Przebieg procesu STFT tak naprawdę nie ma na celu zachowania oryginalnego sygnału.
   • Gdy używasz STFT z nietrywialnym oknem (nie wszystkie), sygnał wejściowy do FFT jest wypaczoną / rozciągniętą wersją oryginalnego fragmentu sygnału.
   • Jest to dobre w przypadku ekstrakcji funkcji, w której niepotrzebne / zbędne dane są odfiltrowywane. Podobnie jak w wykrywaniu sylab, nie wszystkie dane czasowe są wymagane do wykrycia niektórych dźwięków w mowie.
   • Pik w wektorze okna reprezentuje mniejszość pozycji w sygnale audio, na które algorytmy powinny zwrócić uwagę.
2. Zatem surowy wynik odwrotnego STFT może być czymś, czego nie możemy intuicyjnie oczekiwać.
   • Powinny to być okienne fragmenty sygnału, jak powinien wyglądać ifft funkcji STFT.
3. Aby uzyskać oryginalne nieokienne fragmenty sygnału, można zastosować odwrotne okno do surowego wyjścia ifft.
   • Łatwo jest zaprojektować funkcję mapowania, która może cofnąć efekt okna Hanna / Hamminga.
4. Następnie włącza się okno syntezy, aby poradzić sobie z nakładaniem się fragmentacji czasowej
   • Ponieważ oryginalne nieokienne fragmenty sygnału można uznać za już uzyskane, można zastosować dowolne „wagi przejścia” do interpolacji zachodzących na siebie części.
5. Jeśli chcesz wziąć pod uwagę, że fleksja mowy w okienku może mniej szanować słabe sygnały, ale uwielbiać te mocne sygnały, może istnieć sposób zaprojektowania odpowiednich okien syntezy.
6. Można również podać prosty algorytm generowania okna syntezy, stosując następujące zasady:
   • ważą wyższe pozycje w oknie syntezy, jeśli wartość okna analizy dla tej pozycji jest wysoka, w porównaniu z innymi fragmentami pokrywającymi się z tą pozycją.
   • obniżyć wagę pozycji w oknie syntezy, jeśli wartość okna analizy dla tej pozycji jest niska, a inne zachodzące na siebie fragmenty bardziej szanują tę pozycję dzięki większej wartości okna analizy.