Como verificar os resultados da FFT de uma onda senoidal?

Recebi um arquivo de áudio (onda senoidal) 1000Hz como uma entrada para o meu algoritmo FFT. Eu tenho 8192 amostras de espectro de potência em uma matriz.

Qual é a melhor e mais fácil maneira de verificar se minha saída está certa ou errada?

Se eu der um arquivo de áudio silencioso, a saída será zero para todas as amostras. Em uma onda senoidal, o / p aumenta de 20 (0ª amostra) para 26059811 (743ª amostra) e diminui gradualmente para 40.

Se eu tiver uma idéia da faixa de saída, posso provar tecnicamente se a FFT está funcionando.

Qualquer idéia será útil.

Consulte este link para qualquer dúvida técnica.

Parece que você está calculando o espectro calculando a média de 10 janelas (sem sobreposição?) Para obter a magnitude ao quadrado nas frequências 8192 ou 8193 (de 0 a Nyquist, mas alguns algoritmos podem diminuir a frequência de Nyquist na posição 8192).

A primeira coisa a verificar é se o pico está na posição correta. Você não disse qual é a taxa de amostragem, mas o compartimento 743 seria 743/16384 vezes a taxa de amostragem. Se o sinal realmente estiver em 800 Hz, isso coloca os Fs em aproximadamente 17640 amostras / segundo. Isso parece errado. Seu sinal de teste provavelmente estaria em uma taxa padrão como 8000, 16000, 22050, 32000, 44100 ou 48000. Para Fs = 22050, o pico seria muito acentuado na lixeira 800/22050 * 16384 = 594.

Outro critério para verificar é que a energia total no sinal é aproximadamente a mesma nos domínios de tempo e frequência. Aqui está um exemplo em Python:

In [1]: NFFT = 2048; N = 10*NFFT; n = arange(N); Fs = 22050
In [2]: x = 0.4*cos(2*pi*400/Fs*n) + 0.6*cos(2*pi*800/Fs*n)

In [3]: y,freqs = psd(x, NFFT=NFFT, Fs=Fs, pad_to=16384)  # PSD by Welch's Method

In [4]: sum(x**2)/Fs           # time-domain energy
Out[4]: 0.24149869319296949
In [5]: sum(y) * N/16384       # frequency-domain energy
Out[5]: 0.24148752834391252

O sinal de entrada x, que consiste em dois sinusóides amostrados em Fs = 22050 amostras / segundo, é segmentado em 10 janelas não sobrepostas de tamanho NFFT = 2048 amostras. A chamada para psd (densidade espectral de potência) calcula o espectro y como a média da magnitude ao quadrado de dez DFTs de 16384 pontos (na verdade, são 8193 pontos, já que x é valor real).

A energia computada no domínio da frequência tem um fator de escala N / 16384 porque a função psd escalou y para o tamanho da DFT em vez do comprimento total do sinal. Se isso é ou não um problema, depende de como o sistema lida com a normalização do PSD. Outra normalização opcional é a escala de 1 / Fs. Isso combina a energia com o sinal analógico original. As normalizações padrão devem estar bem documentadas na biblioteca.

Você precisa plotar a magnitude da saída da FFT. Eu não estou familiarizado com sua linguagem de programação, mas em Python você usaria algo parecido plot(abs(fft(a))). Para uma entrada silenciosa, a saída deve ser todos os zeros. Para uma entrada de onda senoidal, você verá dois picos:

Para um sinal real, os picos serão simétricos da esquerda para a direita. No entanto, se você estiver realizando uma FFT real (que é mais eficiente em termos de computação), você receberá apenas a metade esquerda do gráfico como sua saída, pois ela ignora a imagem espelhada redundante.

Se a frequência for maior, os picos estarão mais próximos do centro. Se a frequência estiver perfeitamente sincronizada com o tamanho do pedaço , o pico terá apenas um ponto de largura e todo o resto será exatamente 0. Caso contrário, ele terá uma "saia" afilada como acima.

Eu usei a ferramenta de análise de Fourier dentro do Excel Analysis Toolpak para fazer uma verificação rápida dos dados e resultados.