c++ - C ++で周波数スペクトルをプロットする

Question

この質問の下の回答の編集を参照してください。

C++ を使用して正弦波信号の周波数スペクトルをプロットするスクリプトを作成しました。手順は次のとおりです

ハニング窓の適用
fftw3 ライブラリを使用して FFT を適用する

私は3つのグラフを持っています：シグナル、ハニング関数を掛けたときのシグナル、周波数スペクトル。周波数スペクトルが間違っているように見えます。50 Hz にピークがあるはずです。任意の提案をいただければ幸いです。コードは次のとおりです。

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int main()
{
int i;
double y;
int N=50;
double Fs=1000;//sampling frequency
double  T=1/Fs;//sample time 
double f=50;//frequency
double *in;
fftw_complex *out;
double t[N];//time vector 
double ff[N];
fftw_plan plan_forward;

in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

 for (int i=0; i< N;i++)
 {
    t[i]=i*T;
    ff[i]=1/t[i];
    in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
    double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
    in[i] = multiplier * in[i];
  }

  plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );

  fftw_execute ( plan_forward );

  double v[N];

  for (int i = 0; i < N; i++)
    {

    v[i]=20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])/N/2);//Here I have calculated the y axis of the spectrum in dB

    }

   fstream myfile;

   myfile.open("example2.txt",fstream::out);

   myfile << "plot '-' using 1:2" << std::endl;

   for(i = 0; i < N; ++i)

    { 

      myfile << ff[i]<< " " << v[i]<< std::endl;

    }

 myfile.close();

 fftw_destroy_plan ( plan_forward );
 fftw_free ( in );
 fftw_free ( out );
 return 0;
  }

結果をexample2.txtに挿入した後、gnuplotを使用してグラフをプロットしたことを付け加える必要があります。したがって、ff[i] vs v[i] から周波数スペクトルが得られます。

プロットは次のとおりです: それぞれ周波数スペクトルと正弦時間ウィンドウ:

score 1 · Accepted Answer

私の周波数間隔は完全に間違っていました。http://www.ni.com/white-paper/3995/en/#toc1によると; x軸の周波数範囲と分解能は、サンプリングレートとNに依存します。周波数軸の最後の点はFs/2-Fs/Nで、分解能はdF=FS/Nである必要があります。したがって、スクリプトを次のように変更しました。またはサンプリング周波数Fsを下げると、周波数分解能が小さくなり、より良い結果が得られます。)

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int main()
{
int i;
double y;
int N=550;//Number of points acquired inside the window
double Fs=200;//sampling frequency
double dF=Fs/N;
double  T=1/Fs;//sample time 
double f=50;//frequency
double *in;
fftw_complex *out;
double t[N];//time vector 
double ff[N];
fftw_plan plan_forward;

in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

 for (int i=0; i<= N;i++)
 {
 t[i]=i*T;

in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
in[i] = multiplier * in[i];
 }

 for (int i=0; i<= ((N/2)-1);i++)
{ff[i]=Fs*i/N;
}
plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );

fftw_execute ( plan_forward );

double v[N];

for (int i = 0; i<= ((N/2)-1); i++)
{

v[i]=(20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])))/N;  //Here   I  have calculated the y axis of the spectrum in dB

   }

fstream myfile;

myfile.open("example2.txt",fstream::out);

myfile << "plot '-' using 1:2" << std::endl;

for(i = 0;i< ((N/2)-1); i++)

{ 

myfile << ff[i]<< " " << v[i]<< std::endl;

}

 myfile.close();

 fftw_destroy_plan ( plan_forward );
 fftw_free ( in );
 fftw_free ( out );
 return 0;
}

score 0 · Accepted Answer

特に、この Electronics.StackExhcange の投稿を参照してください: https://electronics.stackexchange.com/q/12407/84272 .

50 個のサンプルをサンプリングしているので、25 個の FFT ビンです。1000 Hz でサンプリングしているので、FFT ビンごとに 1000 / 2 / 25 == 250 Hz です。ビンの解像度が低すぎます。

サンプリング周波数を下げるか、サンプル数を増やす必要があると思います。

score -1 · Accepted Answer

あなたの質問がSOにあるので、コードはインデントとスタイルの改善を使用して読みやすくすることができます。

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int main(){
    // use meaningful names for all the variables
    int i;  
    double y;
    int N = 550; // number of points acquired inside the window
    double Fs = 200; // sampling frequency
    double dF = Fs / N;
    double  T = 1 / Fs; // sample time 
    double f = 50; // frequency
    double *in;
    fftw_complex *out;
    double t[N]; // time vector 
    double ff[N];
    fftw_plan plan_forward;

    in = (double*) fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    for (int i = 0; i <= N; i++){
        t[i]=i*T;
        in[i] = 0.7 * sin(2 * M_PI * f * t[i]); // generate sine waveform
        double multiplier = 0.5 * (1 - cos(2 * M_PI * i / (N-1))); // Hanning Window
        in[i] = multiplier * in[i];
    }

    for(int i = 0; i <= ((N/2)-1); i++){
        ff[i] = (Fs * i) / N;
    }

    plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);

    fftw_execute(plan_forward);

    double v[N];
    // Here I have calculated the y axis of the spectrum in dB
    for(int i = 0; i <= ((N/2)-1); i++){
        v[i] = (20 * log(sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]))) / N;  
    }

    fstream myfile;
    myfile.open("example2.txt", fstream::out);
    myfile << "plot '-' using 1:2" << std::endl;

    for(i = 0; i < ((N/2)-1); i++){ 
        myfile << ff[i] << " " << v[i] << std::endl;
    }
    myfile.close();

    fftw_destroy_plan(plan_forward);
    fftw_free(in);
    fftw_free(out);

    return 0;
    }

コードでは、特にループまたは関数呼び出しの前に、より多くのコメントを使用して、入力値 (目的) や戻り値 (結果) を指定できます。

c++ - C ++で周波数スペクトルをプロットする

3 に答える 3

Related

Reference