Interpolação de magnitude da transformada de Fourier discreta (DFT)

7

Por exemplo, para encontrar freqüências de pico, parece válido usar métodos de interpolação com banda limitada nos compartimentos DFT complexos, ou separadamente em suas partes reais e imaginárias, e calcular as magnitudes ou magnitudes ao quadrado dos resultados. Mas e quanto à interpolação limitada por banda das magnitudes dos caixotes do lixo (não acho que isso seja válido) ou de suas magnitudes ao quadrado (talvez válidas)? Por válido, quero dizer que os valores perfeitamente interpolados devem ser iguais aos encontrados, calculando-os a partir de uma DFT maior de uma versão preenchida com zero do sinal no domínio do tempo.

A primeira abordagem garante um resultado não negativo, ao contrário das demais se a interpolação não for perfeita, consulte esta pergunta sobre interpolação não- negativa ou positiva com banda limitada .

Olli Niemitalo
fonte
ei Olli, isso é uma pergunta?
Robert Bristow-johnson
@ robertbristow-johnson Sim, é ... a parte "que tal". Eu tenho minhas suposições lá embora.
Olli Niemitalo 23/10
11
Aliás, o gaussiano é UMA autofunção da Transformada de Fourier. na verdade, há um número infinito de funções próprias e é muito fácil construir uma função própria.
Robert Bristow-johnson
11
@ robertbristow-johnson Ficando um pouco tangencial, mas o que mais é?
Olli Niemitalo 23/10
11
qualquer função de simetria par mais a própria transformação de Fourier (com substituído por ) é uma função própria da Transformada de Fourier. ft
Robert Bristow-johnson

Respostas:

4

Os pontos interpolados da DFT podem ser calculados usando um produto escalar de algumas amostras ao redor da região do pico com um vetor de interpolação pré-calculado. O vetor de interpolação é determinado pela localização da amostra interpolada desejada, levando em consideração a quantidade necessária de preenchimento zero, etc., etc.

Esta técnica e a metodologia para calcular os vetores de interpolação são abordadas no apêndice B deste documento:

http://ericjacobsen.org/FTinterp.pdf

Espero que isto ajude um pouco.

Eric Jacobsen
fonte
Obrigado Eric, é bom ver que você também observou na interpolação de magnitude que "a não linearidade na detecção de magnitude resulta em um sinal que é essencialmente subamostrado". Mostro na minha resposta que a magnitude ao quadrado não é subamostrada quando o sinal no domínio do tempo foi preenchido com zero até 2x. Infelizmente, nem todas as equações no pdf são exibidas corretamente no Chrome, por exemplo, Eq. 3.5
Olli Niemitalo 23/10
11
Ele não é subamostrado até que a magnitude seja detectada, portanto, interpolar os coeficientes complexos para obter o coeficiente interpolado complexo e, em seguida, a magnitude da computação (ou magnitude ao quadrado) evita esse problema.
Eric Jacobsen
Acordado! O pdf funciona bem no Edge.
Olli Niemitalo 23/10
5

Primeiro, uma demonstração de que os quadrados de ambos

[,0 0,0 0,1 1,-1 1,0 0,0 0,] e[,0 0,0 0,1 1,-1 1,0 0,0 0,]

igual

[,0 0,0 0,1 1,-1 1,0 0,0 0,] e

mas os quadrados de suas sinceras interpolações diferem (Fig. 1):

Quadrados de interpolações sinc
Figura 1. Quadrados de interpolações sinc de [1 1,1 1] (azul) e [1 1,-1 1] (vermelho).

Isso demonstra que, em geral, não é possível recuperar o quadrado de um sinal de banda limitada de suas amostras uniformes colhidas na frequência de amostragem crítica do sinal de banda limitada.


Vamos testar diferentes abordagens de interpolação no Octave. O padrão ouro de interpolação com banda limitada é DFT – zero-pad – DFT:

>> format free
>> x = [1 2 3];
>> y = [1 2 3  0 0 0];
>> abs(fft(x))
ans =

 6 1.73205 1.73205

>> abs(fft(y))
ans =

 6 4.3589 1.73205 2 1.73205 4.3589

O último conjunto de números são os valores de magnitude calculados a partir de caixas no domínio da frequência perfeitamente interpoladas. Vamos tentar interpolar a magnitude :

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)))), [0])))
ans =

 6 4.57735 1.73205 0.309401 1.73205 4.57735

Parece um pouco fora. Agora vamos tentar interpolar a magnitude ao quadrado :

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)).^2)), [0])))
ans =

 36 25 3 -8 3 25

>> sqrt(ans)
ans =

 (6,0) (5,0) (1.73205,0) (0,2.82843) (1.73205,0) (5,0)

Não está apenas desativado, mas também sqrt()retornou um número complexo para um valor interpolado negativo. Então, é a única maneira válida de interpolar os valores da posição?

Vamos dar mais uma chance, tentando interpolar os dados do domínio de frequência ampliados por um fator de 2. Por causa da transformação de comprimento uniforme, isso exige a duplicação da "amostra Nyquist", portanto, peça desculpas se o código está ficando difícil de ler.

>> z = [1 2 3  0 0 0  0 0 0  0 0 0];
>> abs(fft(z))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

O acima é o que queremos. Vamos tentar interpolar a magnitude ampliada em 2x :

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)))), [0 0 0]))
ans =

 3.36365 1.10447 0.318175 0 0 0 0 0 0 -0.208949 0.318175 1.10447

>> ans(4) = ans(length(ans)-2)
ans =

 3.36365 1.10447 0.318175 -0.208949 0 0 0 0 0 -0.208949 0.318175 1.10447

>> fft(ans)
ans =

 5.79105 5.59483 4.56785 2.7273 1.5231 1.76882
 2.20895 1.76882 1.5231 2.7273 4.56785 5.59483

Ainda está desligado. Vamos tentar interpolar 2x ao quadrado a amplitude ampliada por amostra :

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)).^2)), [0 0 0]))
ans =

 14 8 3 0 0 0 0 0 0 1.18424e-15 3 8

>> ans(4) = ans(length(ans)-2)
ans =

 14 8 3 1.18424e-15 0 0 0 0 0 1.18424e-15 3 8

>> sqrt(fft(ans))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

Agora funciona perfeito! A mensagem de levar para casa é fazer a ampliação da amostra (zero no domínio do tempo) pelo menos por um fator de dois antes de tentar interpolar no domínio da frequência e interpolar a magnitude ao quadrado em vez da magnitude. Funciona porque assumir a magnitude ao quadrado é o mesmo que multiplicar cada valor de bin pelo seu complexo conjugado. A conjugação complexa preserva a largura de banda da função de banda ilimitada representada pelos dados; portanto, a multiplicação dobra a "largura de banda no domínio do tempo" porque é equivalente à convolução no domínio do tempo. Observe que, ao escolher o método de interpolação, a magnitude quadrada ampliada de 2x ainda é amostrada criticamente; portanto, a amostragem adicional adicional deve facilitar muito a interpolação precisa.

Eu quase respondi minha própria pergunta, mas aceitarei mais informações como resposta!

PS Acabei de descobrir que também existe interpft, o que faz a interpolação com menos sintaxe.

Utilizando informações adicionais

A interpolação torna-se mais fácil ou até exata com informações adicionais sobre os dados, por exemplo, que é um sinc com deslocamento de tempo crítico e amostrado criticamente. Nesse caso, dadas as duas amostrasα pouco antes e γ logo após a maior amostra, o tempo -1 1<d<1 1 do pico pode ser calculado por:

d={0 0E se α=γ,1 1-1 1-p2p,p=γ-αα+γde outra forma,

com o que significa que o sinc é alterado exatamente para o tempo da maior amostra. A mesma interpolação pode ser feita à metade da taxa de amostragem, que é a frequência de amostragem crítica da função subjacente cuja magnitude é igual à de um sinc com desvio no tempo, com as duas amostras sucessivas de maior valor e do quadrado do valor absoluto da função subjacente:d=0 0αβ

0 0d1 1d={0 0E se β=0 0,1 1E se α=0 0,1 1+p-1 1-p22p,p=β-αα+βde outra forma,

As fórmulas não são afetadas pelo dimensionamento da amplitude dos dados. Para a estimativa de frequência, você raramente teria um sinc mudado no tempo puramente real, mas se o fizer, existem fórmulas de interpolação exatas para ele .

Olli Niemitalo
fonte