c++ - C++ 中用于生成频谱图的 vector 和矩阵

Question

这是我第一次尝试使用 C++ 生成正弦信号的频谱图。生成频谱图:

1. 我把真实的正弦信号分成B block

2. 在每个 block 上应用汉宁窗(我假设没有重叠)。这应该给我 fft 的输入, in[j][k]其中 k是区 block 号

3. 申请 fft在 in[j][k]对于每个 block 并存储它。

这是脚本:

#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;
    int N = 500; // sampled
    int Windowsize = 100;
    double Fs = 200; // sampling frequency
    double T = 1 / Fs; // sample time 
    double f = 50; // frequency
    double *in;
    fftw_complex *out;
    double t[N]; // time vector 
    fftw_plan plan_forward;
    std::vector<double> signal(N);
    int B = N / Windowsize; //number of blocks
    in = (double*)fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    //Generating the signal
    for(int i = 0; i < = N; i++){
        t[i] = i * T;
        signal[i] = 0.7 * sin(2 * M_PI * f * t[i]);// generate sine waveform
    }

    //Applying the Hanning window function on each block B
    for(int k = 0; i <= B; k++){ 
        for(int j = 0; j <= Windowsize; j++){
            double multiplier = 0.5 * (1 - cos(2 * M_PI * j / (N-1))); // Hanning Window
            in[j][k] = multiplier * signal[j];
        }
        plan_forward = fftw_plan_dft_r2c_1d (Windowsize, in, out, FFTW_ESTIMATE );
        fftw_execute(plan_forward);
        v[j][k]=(20 * log(sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]))) / N;
    }

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

那么，问题是:声明in[j][k]的正确方法是什么？和 v[j][k]变量。

更新:我已经声明了我的 v [j] [k]作为矩阵:double v [5][249];根据这个网站:http://www.cplusplus.com/doc/tutorial/arrays/所以现在我的脚本看起来像:

#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=500;//Number of pints acquired inside the window
double Fs=200;//sampling frequency
int windowsize=100;
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 tt[5];
double ff[N];
fftw_plan plan_forward;
double  v [5][249];
in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );

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
   }

 for (int k=0; k< 5;k++){ 
   for (int i = 0; i<windowsize; i++){
   double multiplier = 0.5 * (1 - cos(2*M_PI*i/(windowsize-1)));//Hanning Window
   in[i] = multiplier * in[i+k*windowsize];
   fftw_execute ( plan_forward );
         for (int i = 0; i<= (N/2); i++)
         {
         v[k][i]=(20*log10(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i]   [1])));//Here   I  have calculated the y axis of the spectrum in dB
           }
                                    }
                      }

for (int k=0; k< 5;k++)//Center time for each block
      {
       tt[k]=(2*k+1)*T*(windowsize/2);
       }

fstream myfile;
myfile.open("example2.txt",fstream::out);
myfile << "plot '-' using 1:2" << std::endl;
for (int k=0; k< 5;k++){ 
   for (int i = 0; i<= ((N/2)-1); i++)
         { 
        myfile << v[k][i]<< " " << tt[k]<< std::endl;

          }
                         }
myfile.close();
fftw_destroy_plan ( plan_forward );
fftw_free ( in );
fftw_free ( out );
return 0;
  }

我不再收到错误，但频谱图不正确。

Answer 1

如 FFTW's documentation 中所示，使用 fftw_plan_dft_r2c_1d 时输出的大小(在您的情况下为 out)与输入的大小不同。更具体地说，对于 N 个真实样本的输入，输出由 N/2+1 个复数值组成。然后，您可以分配 out :

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));

对于频谱图输出，您将类似地为每个 B block 提供 (N/2+1) 幅度，从而产生二维数组:

double** v = new double*[B];
for (int i = 0; i < B; i++){
  v[i] = new double[(N/2+1)];
}

另请注意，您可以在每次迭代中重用输入缓冲区 in(为新 block 填充数据)。但是，由于您已选择计算 N 点 FFT 并将存储较小的 Windowsize 样本 block (在本例中为 N=500 和Windowsize=100)，确保用零初始化剩余样本:

in = (double*)fftw_malloc(sizeof(double) * N);
for (int i = 0; i < N; i++){
  in[i] = 0;
}

请注意，除了 in 和 v 变量的声明和分配之外，您发布的代码还存在一些其他问题:

• 计算 Hanning window 时，您应该除以 Windowsize-1 而不是 N-1 (因为在您的情况下 N 对应于 FFT 大小)。
• 您正在对同一 block signal 进行一次又一次的 FFT，因为您总是在 [0,Windowsize]< 中使用 j 进行索引范围。您很可能希望在每次处理不同的 block 时都添加一个偏移量。
• 由于 FFT 大小不变，您只需创建一次计划。至少，如果您要在每次迭代中创建计划，您应该在每次迭代中同样销毁它(使用 fftw_destroy_plan)。

还有一些可能需要思考的额外要点:

• 通过除以 N 来缩放对数标度幅度可能不会如您所想。您更有可能想要缩放线性比例幅度(即在取对数之前除以幅度)。请注意，这将导致频谱曲线的恒定偏移，这对于许多应用来说并不那么重要。如果缩放对您的应用程序很重要，您可以查看 another answer of mine了解更多详情。
• 通常用于将线性刻度转换为 decibels 的常用公式 20*log10(x)使用以 10 为底的对数而不是您使用的自然 log (base e~2.7182) 函数。这将导致乘法缩放(拉伸(stretch))，这可能会或可能不会很重要，具体取决于您的应用程序。

总而言之，以下代码可能更符合您的尝试:

// Allocate & initialize buffers
in = (double*)fftw_malloc(sizeof(double) * N);
for (int i = 0; i < N; i++){
  in[i] = 0;
}
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));
v = new (double*)[B];
for (int i = 0; i < B; i++){
  v[i] = new double[(N/2+1)];
}

// Generate the signal
...

// Create the plan once
plan_forward = fftw_plan_dft_r2c_1d (Windowsize, in, out, FFTW_ESTIMATE);

// Applying the Hanning window function on each block B
for(int k = 0; k < B; k++){ 
    for(int j = 0; j < Windowsize; j++){
        // Hanning Window
        double multiplier = 0.5 * (1 - cos(2 * M_PI * j / (Windowsize-1)));
        in[j] = multiplier * signal[j+k*Windowsize];
    }
    fftw_execute(plan_forward);
    for (int j = 0; j <= N/2; j++){
      // Factor of 2 is to account for the fact that we are only getting half
      // the spectrum (the other half is not return by a R2C plan due to symmetry)
      v[k][j] = 2*(out[j][0] * out[j][0] + out[j][1] * out[j][1])/(N*N);
    }
    // DC component and at Nyquist frequency do not have a corresponding symmetric 
    // value, so should not have been doubled up above. Correct those special cases.
    v[k][0]   *= 0.5;
    v[k][N/2] *= 0.5;
    // Convert to decibels
    for (int j = 0; j <= N/2; j++){
      // 20*log10(sqrt(x)) is equivalent to 10*log10(x)
      // also use some small epsilon (e.g. 1e-5) to avoid taking the log of 0
      v[k][j] = 10 * log10(v[k][j] + epsilon);
    }
}

// Clean up
fftw_destroy_plan(plan_forward);
fftw_free(in);
fftw_free(out);
// Delete this last one after you've done something useful with the spectrogram
for (int i = 0; i < B; i++){
  delete[] v[i];
}
delete[] v;