Algoritmo de transformada de Fourier de tiempo corto inverso descrito en palabras

20

Estoy tratando de entender conceptualmente lo que está sucediendo cuando las Transformadas de Fourier de corto tiempo (STFT) de avance e inverso se aplican a una señal de dominio de tiempo discreto. Encontré el artículo clásico de Allen y Rabiner ( 1977 ), así como un artículo de Wikipedia ( enlace ). Creo que también hay otro buen artículo que se puede encontrar aquí .

Estoy interesado en calcular la transformación de Gabor, que no es más que el STFT con una ventana gaussiana.

Esto es lo que entiendo sobre el STFT avanzado :

  1. Se selecciona una subsecuencia de la señal, compuesta de elementos del dominio del tiempo.
  2. La subsecuencia se multiplica por una función de ventana que utiliza la multiplicación punto por punto en el dominio del tiempo.
  3. La subsecuencia multiplicada se toma en el dominio de frecuencia usando la FFT.
  4. Al seleccionar sucesivas subsecuencias superpuestas y repetir el procedimiento anterior, obtenemos una matriz con m filas yn columnas. Cada columna es la subsecuencia calculada en un momento dado. Esto se puede usar para calcular un espectrograma.

Sin embargo, para el STFT inverso , los documentos hablan de una suma sobre las secciones de análisis superpuestas. Me resulta muy difícil visualizar lo que realmente está sucediendo aquí. ¿Qué debo hacer para poder calcular el STFT inverso (en el orden paso a paso como se indicó anteriormente)?

STFT delantero

He creado un dibujo que muestra lo que creo que está sucediendo para el STFT avanzado. Lo que no entiendo es cómo ensamblar cada una de las subsecuencias para que recupere la secuencia de tiempo original. ¿Podría alguien modificar este dibujo o dar una ecuación que muestre cómo se agregan las subsecuencias?Transformar hacia adelante

Transformada inversa

Esto es lo que entiendo sobre la transformación inversa. Cada ventana sucesiva se devuelve al dominio del tiempo utilizando el IFFT. Luego, cada ventana se desplaza por el tamaño del paso y se agrega al resultado del turno anterior. El siguiente diagrama muestra este proceso. La salida sumada es la señal del dominio del tiempo.

Transformación inversa

Ejemplo de código

El siguiente código de Matlab genera una señal de dominio de tiempo sintético, y luego prueba el proceso STFT, demostrando que el inverso es el doble de la transformación directa, dentro del error de redondeo numérico. El principio y el final de la señal están rellenados con ceros para garantizar que el centro de la ventana pueda situarse en el primer y último elemento de la señal en el dominio del tiempo.

N+N01N0

% 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
fuente
2
Gran pregunta, pero, ¿cómo hiciste esos diagramas tan rápido sobre la marcha? ...
Spacey
2
Usé Adobe Illustrator para los diagramas y Mathtype para los caracteres griegos.
Nicholas Kinar
1
"Estoy interesado en calcular la transformación de Gabor, que no es más que el STFT con una ventana gaussiana". Recuerde que la transformación de Gabor es una integral continua, y que las ventanas gaussianas se extienden hasta el infinito. Las implementaciones típicas de STFT usan trozos discretos superpuestos y tienen que usar ventanas de longitud finita.
endolito el
Gracias por señalar eso, endolito. Tiendo a pensar de manera muy discreta cuando hago el procesamiento de la señal.
Nicholas Kinar

Respuestas:

11

El par de transformación STFT se puede caracterizar por 4 parámetros diferentes:

  1. Tamaño FFT (N)
  2. Tamaño de paso (M)
  3. Ventana de análisis (tamaño N)
  4. Ventana de síntesis (tamaño N)

El proceso es el siguiente:

  1. Tome N muestras (tamaño fft) de la ubicación de entrada actual
  2. Aplicar ventana de análisis
  3. Hacer la FFT
  4. Haz lo que quieras hacer en el dominio de frecuencia
  5. FFT inversa
  6. Aplicar ventana de síntesis
  7. Agregar a la salida en la ubicación de salida actual
  8. Ubicación avanzada de entrada y salida por muestras M (tamaño de paso)

El algoritmo de adición de superposición es un buen ejemplo de eso. En este caso, el tamaño del paso es N, el tamaño FFT es 2 * N, la ventana de análisis es rectangular con N unos seguidos de N ceros y las ventanas de síntesis son simplemente todas.

Hay muchas otras opciones para eso y, bajo ciertas condiciones, la transferencia directa / inversa se está reconstruyendo completamente (es decir, puede recuperar la señal original).

El punto clave aquí es que cada muestra de salida generalmente recibe contribuciones aditivas de más de una FFT inversa. La salida debe acumularse en múltiples cuadros. El número de cuadros contribuyentes simplemente viene dado por el tamaño de FFT dividido por el tamaño del paso (redondeado hacia arriba, si es necesario).

Hilmar
fuente
Muchas gracias por su perspicaz respuesta. Entiendo el método de superposición-adición. ¿Qué uso para la ventana de síntesis? ¿Hay una ecuación? Si conozco la función de la ventana de análisis (como una ventana gaussiana), ¿cómo calculo la ventana de síntesis? Entiendo cómo se usa el método de superposición-adición para convolución, pero no entiendo cómo se usa para el STFT. Si el tamaño del paso es step = 1, ¿cómo agrego los marcos juntos? ¿Hay una ecuación?
Nicholas Kinar
Si la función de la ventana de análisis está centrada en cada muestra con un tamaño de paso paso = 1, pongo en cero el principio y el final de la secuencia del dominio del tiempo para que el centro de la ventana esté centrado en cada muestra (incluyendo la primera y la última muestras en la secuencia del dominio del tiempo)?
Nicholas Kinar
Puede elegir ventanas de tamaño de paso, tamaño de fft, análisis y síntesis dependiendo de las necesidades específicas de su aplicación. Un ejemplo es el tamaño de paso N, el tamaño de FFT 2 * N, la evaluación de análisis, la síntesis de todos. Puede modificar eso para análisis sqrt (hanning) y síntesis sqrt (hanning). Cualquiera de los dos funcionará. Me reduzco a lo que haces en el dominio de frecuencia y qué tipo de artefactos como el alias de dominio de tiempo puedes crear.
Hilmar
@Hilmar: Necesito poder hacer modificaciones en el dominio de la frecuencia de la señal, y luego tomar el IFFT para obtener una señal en el dominio del tiempo. Me gustaría minimizar el alias de dominio de tiempo. Todavía no entiendo cómo llevar cada subsecuencia nuevamente al dominio del tiempo y luego sumarlas.
Nicholas Kinar
Escribí un código de prueba y luego actualicé mi pregunta original.
Nicholas Kinar
2

Siete años después de que esta pregunta se planteó por primera vez, me encuentro con esta confusión similar a @Nicholas Kinar. Aquí me gustaría proporcionar algunas ideas y explicaciones perceptivas personales "no oficiales" y "correctas no totalmente aseguradas".

El título de las siguientes declaraciones son exageradas para una mejor inteligibilidad.

  1. El proceso de reenvío de STFT no pretende realmente preservar la señal original.
    • Cuando se usa STFT con una ventana no trivial (No todos), la señal de entrada a FFT es una versión sesgada / estirada del fragmento de señal original.
    • Esto es bueno para la extracción de características, en la que los datos inútiles / redundantes se filtran. Al igual que en la detección de sílabas, no se requieren todos los datos temporales para detectar ciertos tonos en un discurso.
    • El pico en el vector de la ventana representa la minoría de las posiciones en una señal de audio donde los algoritmos deben prestar atención.
  2. Por lo tanto, el resultado bruto del STFT inverso podría ser algo que no podemos esperar intuitivamente.
    • Deben ser los fragmentos de señal en ventana como deberían verse las características de STFT.
  3. Para obtener los fragmentos de señal originales sin ventanas, se puede aplicar una ventana inversa a la salida sin formato de ifft.
    • Es fácil diseñar una función de mapeo que pueda deshacer el efecto de ventana de Hann / Hamming.
  4. La ventana de síntesis se involucra para tratar la superposición de fragmentación temporal
    • Dado que los fragmentos de señal originales sin ventanas pueden verse como ya obtenidos, cualquier "ponderación de transición" puede usarse para interpolar las partes superpuestas.
  5. Si desea tener en cuenta que el fft de un discurso en ventana puede respetar menos las señales débiles pero adorar esas señales potentes, entonces puede haber una manera de diseñar las ventanas de síntesis correspondientes.
  6. Además, se puede proporcionar un algoritmo de generación de ventana de síntesis directo aplicando los siguientes principios:
    • pondera más las posiciones en la ventana de síntesis si el valor de la ventana de análisis para esta posición es alto, en comparación con otros fragmentos que se superponen a esta posición.
    • disminuya las posiciones en la ventana de síntesis si el valor de la ventana de análisis para esta posición es bajo, y otros fragmentos superpuestos respetan más esta posición con un valor de ventana de análisis más grande.
Quirón
fuente
1
Estas son declaraciones interesantes que definitivamente pueden ayudar a fomentar el pensamiento sobre el STFT.
Nicholas Kinar