Algorithme de transformée de Fourier à court terme inverse décrit en mots

20

J'essaie de comprendre conceptuellement ce qui se passe lorsque les transformées de Fourier à court terme avant et inverse (STFT) sont appliquées à un signal discret du domaine temporel. J'ai trouvé l'article classique d'Allen et Rabiner ( 1977 ), ainsi qu'un article Wikipedia ( lien ). Je pense qu'il y a aussi un autre bon article ici .

Je m'intéresse au calcul de la transformée de Gabor, qui n'est rien de plus que la STFT à fenêtre gaussienne.

Voici ce que je comprends du STFT avant :

  1. Une sous-séquence est sélectionnée dans le signal, composée d'éléments du domaine temporel.
  2. La sous-séquence est multipliée par une fonction de fenêtre utilisant une multiplication point par point dans le domaine temporel.
  3. La sous-séquence multipliée est prise dans le domaine fréquentiel à l'aide de la FFT.
  4. En sélectionnant des sous-séquences successives qui se chevauchent et en répétant la procédure ci-dessus, nous obtenons une matrice avec m lignes et n colonnes. Chaque colonne est la sous-séquence calculée à un instant donné. Cela peut être utilisé pour calculer un spectrogramme.

Cependant, pour l' inverse STFT, les articles parlent d'une sommation sur les sections d'analyse qui se chevauchent. Je trouve très difficile de visualiser ce qui se passe vraiment ici. Que dois-je faire pour pouvoir calculer l' inverse STFT (dans l'ordre pas à pas comme ci-dessus)?

Forward STFT

J'ai créé un dessin montrant ce que je pense qui se passe pour le STFT avant. Ce que je ne comprends pas, c'est comment assembler chacune des sous-séquences pour retrouver la séquence temporelle d'origine. Quelqu'un pourrait-il modifier ce dessin ou donner une équation montrant comment les sous-séquences sont ajoutées?Transformation vers l'avant

Transformation inverse

Voici ce que je comprends de la transformation inverse. Chaque fenêtre successive est reprise dans le domaine temporel à l'aide de l'IFFT. Ensuite, chaque fenêtre est décalée de la taille du pas et ajoutée au résultat du décalage précédent. Le diagramme suivant montre ce processus. La sortie sommée est le signal du domaine temporel.

Transformation inverse

Exemple de code

Le code Matlab suivant génère un signal de domaine temporel synthétique, puis teste le processus STFT, démontrant que l'inverse est le double de la transformation directe, avec une erreur d'arrondi numérique. Le début et la fin du signal sont remplis de zéros pour garantir que le centre de la fenêtre peut être situé aux premier et dernier éléments du signal dans le domaine temporel.

N+N0-1N0

% 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);
Nicholas Kinar
la source
2
Grande question - mais, comment avez-vous fait ces diagrammes aussi rapidement à la volée? ...
Spacey
2
J'ai utilisé Adobe Illustrator pour les diagrammes et Mathtype pour les caractères grecs.
Nicholas Kinar le
1
"Je suis intéressé par le calcul de la transformée de Gabor, qui n'est rien de plus que le STFT avec une fenêtre gaussienne." Rappelez-vous que la transformation de Gabor est une intégrale continue et que les fenêtres gaussiennes s'étendent à l'infini. Les implémentations typiques de STFT utilisent des morceaux discrets qui se chevauchent et doivent utiliser des fenêtres de longueur finie.
endolith
Merci de l'avoir signalé, endolith. J'ai tendance à penser de manière très discrète lors du traitement du signal.
Nicholas Kinar

Réponses:

11

La paire de transformations STFT peut être caractérisée par 4 paramètres différents:

  1. Taille FFT (N)
  2. Taille de l'étape (M)
  3. Fenêtre d'analyse (taille N)
  4. Fenêtre de synthèse (taille N)

Le processus est le suivant:

  1. Prenez des échantillons de N (taille fft) à partir de l'emplacement d'entrée actuel
  2. Appliquer la fenêtre d'analyse
  3. Faites la FFT
  4. Faites ce que vous voulez dans le domaine des fréquences
  5. FFT inverse
  6. Appliquer la fenêtre de synthèse
  7. Ajouter à la sortie à l'emplacement de sortie actuel
  8. Avancer l'emplacement d'entrée et de sortie par M (taille de pas) échantillons

L'algorithme d'ajout de chevauchement en est un bon exemple. Dans ce cas, la taille du pas est N, la taille FFT est 2 * N, la fenêtre d'analyse est rectangulaire avec N unités suivies de N zéros et les fenêtres de synthèse sont simplement toutes les unes.

Il y a beaucoup d'autres choix pour cela et dans certaines conditions, le transfert direct / inverse est entièrement reconstruit (c'est-à-dire que vous pouvez récupérer le signal d'origine).

Le point clé ici est que chaque échantillon de sortie reçoit généralement des contributions additives de plus d'une FFT inverse. La sortie doit être accumulée sur plusieurs trames. Le nombre de trames contributrices est simplement donné par la taille de la FFT divisée par la taille du pas (arrondie si nécessaire).

Hilmar
la source
Merci beaucoup pour votre réponse perspicace. Je comprends la méthode de chevauchement-ajout. Qu'est-ce que j'utilise pour la fenêtre de synthèse? Y a-t-il une équation? Si je connais la fonction de fenêtre d'analyse (telle qu'une fenêtre gaussienne), comment puis-je calculer la fenêtre de synthèse? Je comprends comment la méthode de chevauchement-ajout est utilisée pour la convolution, mais je ne comprends pas comment elle est utilisée pour le STFT. Si la taille du pas est step = 1, comment puis-je ajouter les cadres ensemble? Y a-t-il une équation?
Nicholas Kinar le
Si la fonction de fenêtre d'analyse est centrée sur chaque échantillon avec la taille de pas step = 1, dois-je mettre à zéro le début et la fin de la séquence de domaine temporel de sorte que le milieu de la fenêtre soit centré sur chaque échantillon (y compris le premier et le dernier) échantillons dans la séquence du domaine temporel)?
Nicholas Kinar le
Vous pouvez choisir la taille des étapes, la taille des pieds, les fenêtres d'analyse et de synthèse en fonction des besoins spécifiques de votre application. Un exemple est la taille de pas N, la taille FFT 2 * N, la suspension d'analyse, la synthèse de tous. Vous pouvez modifier cela en sqrt d'analyse (suspension) et sqrt de synthèse (suspension). L'un ou l'autre fonctionnera. Je me résume à ce que vous faites dans le domaine fréquentiel et au type d'artefacts tels que l'alias de domaine temporel que vous pouvez créer.
Hilmar
@ Hilmar: J'ai besoin de pouvoir apporter des modifications dans le domaine fréquentiel au signal, puis prendre l'IFFT pour obtenir un signal dans le domaine temporel. Je voudrais minimiser l'alias de domaine temporel. Je ne comprends toujours pas comment reprendre chaque sous-séquence dans le domaine temporel puis les ajouter ensemble.
Nicholas Kinar le
J'ai écrit du code de test, puis mis à jour ma question d'origine.
Nicholas Kinar le
2

Sept ans après que cette question a été soulevée pour la première fois, je rencontre cette confusion similaire à @Nicholas Kinar. Ici, je voudrais fournir quelques idées et explications personnelles «non officielles» et «exactitude pas entièrement assurées».

Le titre des déclarations suivantes est exagéré pour une meilleure intelligibilité.

  1. Le processus direct de STFT n'est pas vraiment destiné à conserver le signal d'origine.
    • Lorsque vous utilisez STFT avec une fenêtre non triviale (pas tout-en-un), le signal d'entrée vers FFT est une version asymétrique / étirée du fragment de signal d'origine.
    • C'est bon pour l'extraction de fonctionnalités, dans laquelle les données inutiles / redondantes sont filtrées. Comme dans la détection des syllabes, toutes les données temporelles ne sont pas nécessaires pour détecter certaines tonalités dans un discours.
    • Le pic du vecteur fenêtre représente la minorité des positions d'un signal audio auxquelles les algorithmes doivent prêter attention.
  2. Ainsi, le résultat brut de l'inverse STFT pourrait être quelque chose que nous ne pouvons pas attendre intuitivement.
    • Ce devrait être les fragments de signaux fenêtrés à quoi devraient ressembler les différentes fonctionnalités de STFT.
  3. Pour obtenir les fragments de signal sans fenêtre d'origine, on peut appliquer une fenêtre inverse à la sortie brute de ifft.
    • Il est facile de concevoir une fonction de cartographie qui peut annuler l'effet de fenêtre hann / hamming.
  4. La fenêtre de synthèse est alors impliquée pour faire face au chevauchement de la fragmentation temporelle
    • Etant donné que les fragments de signal non fenêtrés d'origine peuvent être considérés comme déjà obtenus, toute "pondération de transition" peut être utilisée pour interpoler les parties superposées.
  5. Si vous souhaitez considérer que la fft d'un discours fenêtré peut moins respecter les signaux faibles mais adorer ces signaux puissants, alors il peut y avoir un moyen de concevoir des fenêtres de synthèse correspondantes.
  6. En outre, un algorithme de génération de fenêtre de synthèse simple peut être donné en appliquant les principes suivants:
    • pondérer plus les positions dans la fenêtre de synthèse si la valeur de la fenêtre d'analyse pour cette position est élevée, en comparaison avec d'autres fragments qui chevauchent cette position.
    • poids abaisser les positions dans la fenêtre de synthèse si la valeur de la fenêtre d'analyse pour cette position est faible, et d'autres fragments qui se chevauchent honorent davantage cette position avec une valeur de fenêtre d'analyse plus grande.
Chiron
la source
1
Ce sont des déclarations intéressantes qui peuvent certainement aider à encourager la réflexion sur le STFT.
Nicholas Kinar