math - 正弦波周波数フィッティング

Question

この質問は、以前の同様の質問に基づいています。

次の式と調整済み (いくつかのランダム データ) があります: 0.44*sin(N* 2*PI/30)

FFT を使用して、生成されたデータから周波数を取得しようとしています。ただし、周波数は最終的に周波数に近くなりますが、等しくはなりません (これにより、意図したよりも波が少し大きくなります)。

FFT の最大周波数は 7hz ですが、予想される周波数は (30/2PI) 4.77hz です。

FFT とプロットされた値のグラフを含めました。

私が使用しているコードは次のとおりです。

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

正の FFT はここにあります。基本的に、FFT グラフを中央に配置し、負の信号をカットします。

私の質問は、周波数だけの最小二乗法に頼らずに FFT をより正確にするにはどうすればよいかということです。

Answer 1

FFTは、（準）周期信号の高解像度周波数測定には適していないと思います。以下を参照してください。

すべての離散FFTは、非整数のビン周波数（つまり、特定のFFTの周波数ステップの1つに正確に対応しない任意の周波数）に広がります。これらの「中間」周波数は、最も近い整数ビンの周りに塗り付け/拡散されます。この広がりの形状（「広がり関数」）は、FFTに使用されるウィンドウ関数によって異なります。この拡散関数（物事を単純化および一般化するため）は、非常に狭いが非常に不規則（非常に高い山/非常に低い谷）であるか、または広くても不規則ではありません。理論的には、正弦波の非常に細かい周波数掃引を実行し、それぞれのFFTを計算できます。次に、すべてのFFTの出力を、その出力を生成した周波数と一緒に保存することで、関数の形状と動作を「較正」できます。

たくさんの努力。

ただし、単一の信号の周波数のみを測定する必要がある場合は、これを行わないでください。

代わりに、波長を測定してみてください。これは、サンプルのゼロクロス間の距離を測定し（おそらく、より正確にするために複数のサイクルの場合、1000サイクルを測定します）、サンプルレートをそれで割って周波数に到達するのと同じくらい簡単です。はるかに単純で、高速で、はるかに正確です。

例：48000 Hzのサンプルレート、4.77 Hzの信号は、最も粗いアプローチで1サイクルの長さを測定するだけで、約0.0005Hzの分解能になります。（nサイクルかかる場合、周波数分解能もn倍になります。）

Answer 2

他の人が述べたように、信号の周波数を誤解しています。いくつかのことを明確にするために例を挙げましょう。

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')

これで、FFT のピークが 6Hz の信号の元の周波数に対応することがわかります。

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

Parseval の定理が成り立つことも確認できます。

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

---

オプション 2

または、信号処理ツールボックスの関数を使用できます。

%# 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)

この関数は対数スケールを使用することに注意してくださいplot(fr,10*log10(Pxx)):plot(fr,Pxx)

Answer 3

あなたが探しているのは周波数推定法で、たくさんあります。FFT は、いくつかの推定方法の 1 つのコンポーネントです。あなたの例のように、ピーク振幅ビンを使用するだけで、最悪の解像度が得られます（ただし、他の正確に周期的な正弦波に対するノイズ耐性は最大になります）。ノイズの少ない状況では、補間できます。対数振幅の放物線補間は一般的な推定量の 1 つですが、矩形ウィンドウの場合は FFT 結果の同期補間の方が適している場合があります。ゼロパディングとより長いFFTの実行は、基本的に補間と同等です。

ゼロ ノイズで正確な正弦波を得るには、FFT を忘れて、3 つの未知数で方程式を解くだけです。これには、エイリアス化されていないサンプル ポイントが 3 つまたは 4 つ含まれる場合があります。

私のDSP Web ページに、他のいくつかの周波数推定方法をリストしています。

Answer 4

N を秒単位の時間とすると、周波数は 1/30Hz ( y=A * sin( 2* PI * f * t))

周波数分解能 = サンプルレート / FFT ポイント

サンプル レートはナイキスト基準によって決定されます。サンプル レート (サンプル/秒) は、分析する最大周波数の少なくとも 2 倍である必要があります。たとえば、24kHz までの分析では 48kHz です。(「実際の」データの場合は、少しバッファを用意しておくとよいでしょう)。

そのため、FFT のサイズを大きくする必要がある場合があります。

Answer 5

関数から生成する場合は、サンプルを操作するのではなく、多くのポイントを生成してBIG fftを実行できるため、周波数ビンは非常に小さく、高精度になります。しかし、それは基本的な問題を解決しません。

Answer 6

まず、質問の訂正: (30/2PI) は周波数ではありません。信号の周波数は、使用したサンプリング レートの 1/30 * です。次に、sampledata ベクトルの長さを教えてください。FFT が値のベクトルを返す場合、i 番目の値は f_i = i/N に対応します。ここで、N はベクトルの長さで、i \in [0,N-1]私。つまり、N は 30*i に等しくなければなりません。つまり、N は 30 の倍数でなければなりません。使用したベクトルの長さは 30 の倍数でしたか? そうでない場合は、それを作成してみてください。それで問題は解決するはずです。

Answer 7

ウィンドウ機能を試してみませんか？