Algoritmo inverso de transformada de Fourier de tempo curto descrito em palavras

20

Estou tentando entender conceitualmente o que está acontecendo quando as Transformadas de Fourier de Tempo Curto (STFT), inversas e inversas, são aplicadas a um sinal discreto no domínio do tempo. Encontrei o artigo clássico de Allen e Rabiner ( 1977 ), bem como um artigo da Wikipedia ( link ). Acredito que também haja outro bom artigo aqui .

Estou interessado em calcular a transformação Gabor, que nada mais é do que o STFT com uma janela gaussiana.

Isto é o que eu entendo sobre o STFT avançado :

  1. Uma sub-sequência é selecionada a partir do sinal, composto por elementos no domínio do tempo.
  2. A subseqüência é multiplicada por uma função de janela usando multiplicação ponto a ponto no domínio do tempo.
  3. A sub-sequência multiplicada é inserida no domínio da frequência usando a FFT.
  4. Selecionando sucessivas sub sequências sobrepostas e repetindo o procedimento acima, obtemos uma matriz com m linhas e n colunas. Cada coluna é a sub-sequência calculada em um determinado momento. Isso pode ser usado para calcular um espectrograma.

No entanto, para o STFT inverso , os artigos falam sobre um somatório sobre as seções de análise sobrepostas. Estou achando muito desafiador visualizar o que realmente está acontecendo aqui. O que tenho que fazer para poder calcular o STFT inverso (em ordem passo a passo, como acima)?

Encaminhar STFT

Eu criei um desenho mostrando o que acho que está acontecendo no STFT avançado. O que não entendo é como montar cada uma das sub-sequências para que eu volte à sequência de tempo original. Alguém poderia modificar esse desenho ou fornecer uma equação mostrando como as sub sequências são adicionadas?Transformação direta

Transformação Inversa

Aqui está o que eu entendo sobre a transformação inversa. Cada janela sucessiva é retornada ao domínio do tempo usando o IFFT. Em seguida, cada janela é alterada pelo tamanho da etapa e adicionada ao resultado da mudança anterior. O diagrama a seguir mostra esse processo. A saída somada é o sinal no domínio do tempo.

Transformação inversa

Exemplo de código

O código Matlab a seguir gera um sinal sintético no domínio do tempo e, em seguida, testa o processo STFT, demonstrando que o inverso é o duplo da conversão direta, com erro de arredondamento numérico. O início e o fim do sinal são preenchidos com zero para garantir que o centro da janela possa estar situado no primeiro e no último elementos do sinal no domínio do tempo.

N+N0 0-1N0 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);
Nicholas Kinar
fonte
2
Grande questão -, mas, como você fazer esses diagramas isso rapidamente na mosca ...?
Spacey
2
Usei o Adobe Illustrator para os diagramas e Mathtype para os caracteres gregos.
Nicholas Kinar
1
"Estou interessado em calcular a transformação Gabor, que nada mais é do que o STFT com uma janela gaussiana". Lembre-se de que a transformação de Gabor é uma integral contínua e que as janelas gaussianas se estendem ao infinito. As implementações típicas do STFT usam pedaços discretos sobrepostos e precisam usar janelas de comprimento finito.
Endolith
Obrigado por apontar isso, endólito. Costumo pensar de uma maneira muito discreta ao processar o sinal.
Nicholas Kinar

Respostas:

11

O par de transformações STFT pode ser caracterizado por 4 parâmetros diferentes:

  1. Tamanho FFT (N)
  2. Tamanho da etapa (M)
  3. Janela de análise (tamanho N)
  4. Janela de síntese (tamanho N)

O processo é como se segue:

  1. Pegue amostras N (tamanho de pés) do local de entrada atual
  2. Janela de análise de aplicação
  3. Faça a FFT
  4. Faça o que quiser no domínio da frequência
  5. FFT inversa
  6. Aplicar janela de síntese
  7. Adicione à saída no local de saída atual
  8. Localização avançada de entrada e saída por amostras M (tamanho da etapa)

O algoritmo de sobreposição de adição é um bom exemplo disso. Nesse caso, o tamanho da etapa é N, o tamanho da FFT é 2 * N, a janela de análise é retangular com N, seguida de N zeros e as janelas de síntese são simplesmente todas.

Existem muitas outras opções para isso e, sob certas condições, a transferência direta / inversa está sendo totalmente reconstruída (ou seja, você pode recuperar o sinal original).

O ponto principal aqui é que cada amostra de saída geralmente recebe contribuições aditivas de mais de uma FFT inversa. A saída precisa ser acumulada em vários quadros. O número de quadros contribuintes é simplesmente dado pelo tamanho da FFT dividido pelo tamanho da etapa (arredondado para cima, se necessário).

Hilmar
fonte
Muito obrigado pela sua resposta perspicaz. Eu entendo o método de sobreposição de adição. O que eu uso para a janela de síntese? Existe uma equação? Se conheço a função da janela de análise (como uma janela gaussiana), como computo a janela de síntese? Entendo como o método de sobreposição-adição é usado para convolução, mas não entendo como é usado para o STFT. Se o tamanho da etapa for etapa = 1, como adiciono os quadros? Existe uma equação?
Nicholas Kinar
Se a função da janela de análise estiver centralizada em cada amostra com o tamanho da etapa step = 1, zero o início e o final da sequência no domínio do tempo para que o meio da janela seja centralizado em cada amostra (incluindo a primeira e a última) amostras na sequência no domínio do tempo)?
Nicholas Kinar
Você pode escolher o tamanho da etapa, tamanho do fft, janelas de análise e síntese, dependendo das necessidades específicas do seu aplicativo. Um exemplo é o tamanho da etapa N, tamanho da FFT 2 * N, análise, execução, síntese. Você pode modificar isso para análise sqrt (hanning) e síntese sqrt (hanning). Qualquer um funcionará. Eu me refiro ao que você faz no domínio da frequência e a que tipo de artefatos, como alias do domínio do tempo, você pode criar.
Hilmar
@ Hilmar: Eu preciso fazer modificações no domínio da frequência no sinal e, em seguida, pegar o IFFT para obter um sinal no domínio do tempo. Gostaria de minimizar o alias do domínio do tempo. Ainda não entendo como levar cada subseqüência de volta ao domínio do tempo e depois adicioná-las.
Nicholas Kinar
Escrevi um código de teste e atualizei minha pergunta original.
Nicholas Kinar
2

Sete anos depois que essa pergunta foi levantada, eu me deparei com essa confusão semelhante a @Nicholas Kinar. Aqui eu gostaria de fornecer algumas idéias e explicações perceptivas pessoais "não oficiais" e "corretivas não totalmente garantidas".

O título das declarações a seguir é exagerado para melhor inteligibilidade.

  1. O processo de encaminhamento do STFT não tem como objetivo preservar o sinal original.
    • Ao usar STFT com uma janela não trivial (nem todas), o sinal de entrada para FFT é uma versão distorcida / esticada do fragmento de sinal original.
    • Isso é bom para a extração de recursos, na qual os dados inúteis / redundantes são filtrados. Como na detecção de sílabas, nem todos os dados temporais são necessários para detectar certos tons em um discurso.
    • O pico no vetor da janela representa a minoria das posições em um sinal de áudio em que os algoritmos devem prestar atenção.
  2. Portanto, o resultado bruto do STFT inverso pode ser algo que talvez não esperemos intuitivamente.
    • Devem ser os fragmentos de sinal em janela com os quais o ifft dos recursos STFT deve se parecer.
  3. Para obter os fragmentos originais de sinal sem janela, pode-se aplicar uma janela inversa à saída bruta do ifft.
    • É fácil projetar uma função de mapeamento que pode desfazer o efeito da janela hann / hamming.
  4. A janela de síntese é então envolvida para lidar com a sobreposição da fragmentação temporal
    • Como os fragmentos originais de sinal sem janela podem ser vistos como já obtidos, quaisquer "ponderações de transição" podem ser usadas para interpolar as partes sobrepostas.
  5. Se você gostaria de considerar que o tom de uma fala em janela pode respeitar menos os sinais fracos, mas adora esses sinais poderosos, pode haver uma maneira de projetar as janelas de síntese correspondentes.
  6. Além disso, um algoritmo de geração de janela de síntese direta pode ser fornecido aplicando os seguintes princípios:
    • pesar mais alto as posições na janela de síntese se o valor da janela de análise para essa posição for alto, comparado com outros fragmentos que se sobrepõem a essa posição.
    • o peso diminui as posições na janela de síntese se o valor da janela de análise para essa posição for baixo e outros fragmentos sobrepostos honram mais essa posição com um valor maior da janela de análise.
Quíron
fonte
1
Estas são declarações interessantes que podem definitivamente ajudar a incentivar o pensamento sobre o STFT.
Nicholas Kinar 17/11