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 :
- Se selecciona una subsecuencia de la señal, compuesta de elementos del dominio del tiempo.
- La subsecuencia se multiplica por una función de ventana que utiliza la multiplicación punto por punto en el dominio del tiempo.
- La subsecuencia multiplicada se toma en el dominio de frecuencia usando la FFT.
- 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?
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.
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.
% 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);
fuente
Respuestas:
El par de transformación STFT se puede caracterizar por 4 parámetros diferentes:
El proceso es el siguiente:
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).
fuente
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.
fuente