Ajuste de frequência de onda seno

Question

Esta pergunta é baseada em um Pergunta semelhante anterior.

Eu tenho a seguinte equação e um ajustado (alguns dados aleatórios): 0,44*sin (n*2*pi/30)

Estou tentando usar o FFT para obter a frequência dos dados gerados. No entanto, a frequência acaba sendo próxima, mas não igual à frequência (o que torna a onda um pouco maior do que o pretendido)

As frequências que estão no máximo para a FFT são 7Hz, no entanto, a frequência esperada é (30/2pi) 4,77Hz.

Incluí um gráfico dos valores FFT e plotados.

O código que estou usando é:

[sampleFFTValues sFreq] = positiveFFT(sampledata, 1);
sampleFFTValues = abs(sampleFFTValues);
[v sFFTV]= max(sampleFFTValues)

FFT positivo pode ser encontrado aqui. Basicamente, ele centra o gráfico da FFT e corta os sinais negativos.

Minha pergunta é como posso fazer com que a FFT seja mais precisa sem ter que recorrer a mínimos quadrados apenas para a frequência?

Answer 1

Não acho que a FFT seja boa para uma medição de frequência de resolução fina para sinais periódicos (quase) - veja abaixo.

Todo FFT discreto tem se espalhando em frequências de bandas não inteiras (que está em qualquer frequência que não corresponda exatamente a uma das etapas de frequência da FFT específica); Essas frequências "intermediárias" serão manchadas/espalhadas pela lixeira inteira mais próxima. A forma dessa propagação ("Função de espalhamento") depende da função de janela usada para a FFT. Essa função de espalhamento - para simplificar e generalizar as coisas - é muito estreita, mas muito esfarrapada (picos muito altos/vales muito baixos) ou mais largos, mas menos esfarrapados. Em teoria, você pode fazer uma varredura de frequência muito fina de ondas senoidais e calcular a FFT para cada uma delas, e então você pode "calibrar" a forma e o comportamento da função salvando savendo saveiros de todas as FFTs juntamente com a frequência que resultou nessa saída, e depois comparando a saída FFT do sinal a ser medido com os resultados salvos anteriormente e encontrando o "mais próximo", encontre uma frequência mais exata.

Muito esforço.

Mas não faça isso se você precisar apenas medir a frequência de um único sinal.

Em vez disso, tente medir o comprimento de onda. Isso pode ser tão simples quanto medir a distância entre zero cruzamentos (talvez para vários ciclos obter mais precisão - diabos, meça 1000 ciclos, se você tiver tantos) em amostras e divida a taxa de amostragem por chegar à frequência. Muito mais simples, mais rápido e muito mais preciso.

Exemplo: taxa de amostragem de 48000 Hz, sinal de 4,77 Hz resulta em resolução de ~ 0,0005 Hz apenas medindo o comprimento de 1 Ciclo com a abordagem mais cruel. (Se você pegar n ciclos, a resolução de frequência se multiplica por n também.)

Answer 2

Como mencionado por outros, você está interpretando mal a frequência do sinal. Deixe -me dar um exemplo para esclarecer algumas coisas:

Fs = 200;                        %# sampling rate
t = 0:1/Fs:1-1/Fs;               %# time vector of 1 second 
f = 6;                           %# frequency of signal
x = 0.44*sin(2*pi*f*t);          %# sine wave

N = length(x);                   %# length of signal
nfft = N;                        %# n-point DFT, by default nfft=length(x)
                                 %# (Note: it is faster if nfft is a power of 2)
X = abs(fft(x,nfft)).^2 / nfft;  %# square of the magnitude of FFT

cutOff = ceil((nfft+1)/2);       %# nyquist frequency
X = X(1:cutOff);                 %# FFT is symmetric, take first half
X(2:end -1) = 2 * X(2:end -1);   %# compensate for the energy of the other half
fr = (0:cutOff-1)*Fs/nfft;       %# frequency vector

subplot(211), plot(t, x)
title('Signal (Time Domain)')
xlabel('Time (sec)'), ylabel('Amplitude')

subplot(212), stem(fr, X)
title('Power Spectrum (Frequency Domain)')
xlabel('Frequency (Hz)'), ylabel('Power')

Agora você pode ver que o pico na FFT corresponde à frequência original do sinal em 6Hz

[v idx] = max(X);
fr(idx)
ans = 
      6

Podemos até verificar isso Teorema de Parseval Holds:

( sum(x.^2) - sum(X) )/nfft < 1e-6

---

opção 2

Como alternativa, podemos usar as funções da caixa de ferramentas de processamento de sinal:

%# estimate the power spectral density (PSD) using the periodogram
h = spectrum.periodogram;
hopts = psdopts(h);
set(hopts, 'Fs',Fs, 'NFFT',nfft, 'SpectrumType','onesided')

hpsd = psd(h, x, hopts);
figure, plot(hpsd)

Pxx = hpsd.Data;
fr = hpsd.Frequencies;
[v idx]= max(Pxx)
fr(idx)

avgpower(hpsd)

Observe que esta função usa uma escala logarítmica: plot(fr,10*log10(Pxx)) ao invés de plot(fr,Pxx)

Answer 3

Assumindo que n é tempo em segundos, sua frequência é de 1/30Hz (y=A * sin( 2* PI * f * t))

Resolução de frequência = taxa de amostra / ponto de FFT

A taxa de amostragem é determinada pelo critério nyquist, a taxa de amostra (amostras/segundo) deve ser pelo menos duas vezes a frequência máxima a ser analisada, por exemplo, 48kHz para analisar até 24kHz. (Para dados da "vida real", é bom ter um buffer um pouco).

Portanto, pode ser necessário aumentar o tamanho da sua FFT.

Answer 4

O que você está procurando é um método de estimativa de frequência e há muitos. Um FFT é um componente de vários métodos de estimativa. Apenas usando o compartimento de magnitude de pico, como no seu exemplo, oferece a pior resolução (mas a maior imunidade de ruído para qualquer outro sinusóide exatamente periódico). Em situações de baixo ruído, você pode interpolar. A interpolação parabólica da magnitude do log é um estimador comum, mas a interpolação de sincronização dos resultados da FFT pode ser melhor para uma janela retangular. O acalmar zero e fazer uma FFT mais longa é basicamente equivalente à interpolação.

Para um sinusóide exato em ruído zero, esqueça a FFT e basta resolver a equação em 3 incógnitas, o que pode envolver apenas 3 ou 4 pontos de amostra não perseguidos, algoritmos para fazer isso aqui e aqui.

Eu listo alguns outros métodos de estimativa de frequência no meu Página da Web do DSP.

Answer 5

Se você estiver gerando a partir de uma função, em vez de trabalhar com amostras, poderá gerar muitos pontos e executar uma grande FFT para que as caixas de frequência sejam muito pequenas para alta precisão. Mas não resolverá o problema básico.

Answer 6

Primeiro, uma correção para sua pergunta: (30/2PI) não é a frequência. A frequência do seu sinal é 1/30 * qualquer taxa de amostragem que você usou. Segundo, você pode me dizer qual era o comprimento do vetor de sampledata? Quando a FFT retorna um vetor de valores, o valor é eu. Em outras palavras, n deve ser igual a 30*i, ou seja, n deve ser um múltiplo de 30. Agora, o comprimento do vetor que você usou, um múltiplo de 30? Se não tentar fazê -lo, e isso deve resolver o problema.

Answer 7

Tente um função de janela?