FFT結果の大きさは波の周波数に依存しますか?
-
21-09-2019 - |
質問
私はFFTから得た結果に困惑しており、どんな助けにも感謝しています。
FFTW 3.2.2を使用していますが、他のFFT実装(Java)と同様の結果が得られています。正弦波のFFTを採取すると、結果のスケーリングは、波の周波数(Hz)に依存します。特に、整数に近いかどうかにかかわらず。結果の値は、周波数が整数に近い場合に非常に小さく、頻度が整数の間にあるときに数桁大きくなります。 このグラフ さまざまな周波数の波の周波数に対応するFFT結果のスパイクの大きさを示します。これは正しいですか??
FFTの逆FFTは、サンプルの数が元の正弦波時間に等しいことを確認しました。 FFTの形状も正しいようです。
個々の正弦波を分析していても、それほど悪くはないでしょう。なぜなら、その高さに関係なくFFTのスパイクを探すことができるからです。問題は、正弦波の合計を分析したいということです。たとえば、440 Hzと523.25 Hzで正弦波の合計を分析すると、523.25 Hzのスパイクのみが表示されます。もう一方のスパイクは非常に小さく、ノイズのように見えます。 Matlabでは機能しているため、この作業を行うには何らかの方法が必要です。両方の周波数で同様のサイズのスパイクを取得します。以下のコードを変更して、さまざまな周波数のスケーリングを均等化するにはどうすればよいですか?
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <fftw3.h>
#include <cstdio>
using namespace std;
const double PI = 3.141592;
/* Samples from 1-second sine wave with given frequency (Hz) */
void sineWave(double a[], double frequency, int samplesPerSecond, double ampFactor);
int main(int argc, char** argv) {
/* Args: frequency (Hz), samplesPerSecond, ampFactor */
if (argc != 4) return -1;
double frequency = atof(argv[1]);
int samplesPerSecond = atoi(argv[2]);
double ampFactor = atof(argv[3]);
/* Init FFT input and output arrays. */
double * wave = new double[samplesPerSecond];
sineWave(wave, frequency, samplesPerSecond, ampFactor);
double * fftHalfComplex = new double[samplesPerSecond];
int fftLen = samplesPerSecond/2 + 1;
double * fft = new double[fftLen];
double * ifft = new double[samplesPerSecond];
/* Do the FFT. */
fftw_plan plan = fftw_plan_r2r_1d(samplesPerSecond, wave, fftHalfComplex, FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(plan);
memcpy(fft, fftHalfComplex, sizeof(double) * fftLen);
fftw_destroy_plan(plan);
/* Do the IFFT. */
fftw_plan iplan = fftw_plan_r2r_1d(samplesPerSecond, fftHalfComplex, ifft, FFTW_HC2R, FFTW_ESTIMATE);
fftw_execute(iplan);
fftw_destroy_plan(iplan);
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < samplesPerSecond; i++) {
printf("\t%.6f", wave[i]);
}
printf("\n");
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < fftLen; i++) {
printf("\t%.9f", fft[i]);
}
printf("\n");
printf("\n");
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < samplesPerSecond; i++) {
printf("\t%.6f (%.6f)", ifft[i], samplesPerSecond * wave[i]); // actual and expected result
}
delete[] wave;
delete[] fftHalfComplex;
delete[] fft;
delete[] ifft;
}
void sineWave(double a[], double frequency, int samplesPerSecond, double ampFactor) {
for (int i = 0; i < samplesPerSecond; i++) {
double time = i / (double) samplesPerSecond;
a[i] = ampFactor * sin(2 * PI * frequency * time);
}
}
解決
結果の値は、周波数が整数に近い場合に非常に小さく、頻度が整数の間にあるときに数桁大きくなります。
これは、高速フーリエ変換が入力が周期的であり、無限に繰り返されると想定しているためです。サイン波の非統合数があり、この波形を繰り返す場合、それは完全な正弦波ではありません。これにより、fft結果が生じます 「スペクトル漏れ」
調べてください ウィンドウ関数. 。これらは、最初と最後に入力を減衰させるため、スペクトル漏れが減少します。
PS:基本的に正確な周波数コンテンツを取得したい場合は、多くの波のサイクルをキャプチャし、サイクルあたりあまり多くのポイントをキャプチャする必要はありません(サイクルあたり32ポイントまたは64ポイントはおそらく十分です)。より高い高調波で正確な周波数コンテンツを取得する場合は、サイクルごとに少数のサイクル数をキャプチャします。
他のヒント
GNUラジオコードを見ることをお勧めします。あなたにとって特に興味深いファイルはusrp_fft.pyです。