6

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); 
 }
}
4

2 に答える 2

10

結果の値は、周波数が整数に近い場合は非常に小さくスケーリングされ、周波​​数が整数の間にある場合は桁違いに大きくなります。

これは、高速フーリエ変換が入力が周期的であり、無限に繰り返されると想定しているためです。正弦波の数が整数でなく、この波形を繰り返す場合、それは完全な正弦波ではありません。これにより、FFT の結果に「スペクトル漏れ」が発生します。

ウィンドウ関数を調べます。これらは、最初と最後で入力を減衰させるため、スペクトルの漏れが減少します。

追記: 基本波周辺の正確な周波数コンテンツを取得したい場合は、多くの波のサイクルをキャプチャし、サイクルごとにあまり多くのポイントをキャプチャする必要はありません (サイクルごとに 32 または 64 ポイントでおそらく十分です)。高調波で正確な周波数コンテンツを取得したい場合は、キャプチャするサイクル数を減らし、サイクルごとにより多くのポイントを取得します。

于 2009-12-24T22:16:04.160 に答える
0

GNU Radio コードを確認することをお勧めします。特に興味のあるファイルは usrp_fft.py です。

于 2009-12-24T21:45:02.197 に答える