Ajuste de frequência de onda seno
-
22-09-2019 - |
Pergunta
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?
Solução
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.)
Outras dicas
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)
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.
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.
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.
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.
Tente um função de janela?