c - 需要帮助使这个程序计算 PSD

Question

我正在尝试制作一个程序来计算时间序列的 PSD(16384 个样本)，假设这里的正弦曲线是代码:

// generating sin samples 

#include <stdio.h>
#include <math.h>
#include <complex.h>

int main(){
    FILE *inp =NULL,*inp2;
    double value = 0.0; 
    float frequency = 1; // signal frequency 
    double timeSec = 1 ; // time in sec 
    unsigned int  numberOFSamples = 2048*8;
    double steps = timeSec / numberOFSamples; 
    double timer =0.0;
    float  dcValue =0.0;
    double index = 0;
     inp = fopen("sinus","wb+");
     inp2= fopen("sinusD","wb+");
    for( timer=0.0 ; timer<timeSec;timer+=steps){
        value= sin(2*M_PI*frequency*timer) +dcValue;
        fprintf(inp,"%lf ",value);
        fwrite(&value,sizeof(double),1,inp2);

     }

     fclose(inp);
     fclose(inp2);
     return 0;

}

正弦值的生成是正确的我现在使用 Matlab 检查它 PSD 应该是 1024 大这里是代码:

    #include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <complex.h>

#define WINDOW_SIZE 1024

int main (){
  FILE* inputFile = NULL;
  FILE* outputFile= NULL;
  double* inputData=NULL; 
  double* outputData=NULL;
  double* windowData=NULL;
  unsigned int windowSize = 1024;
  int overlaping =0;
  int index1 =0,index2=0, i=0;
  double powVal= 0.0;
  fftw_plan plan_r2hc;
  double w[WINDOW_SIZE];
  for (i=0; i<WINDOW_SIZE; i++) {
  w[i] = (1.0 - cos(2.0 * M_PI * i/(WINDOW_SIZE-1))) * 0.5;  // hann window
  }

// mememory allocation
  inputData = (double*) fftw_malloc(sizeof(double)*windowSize);
  outputData= (double*) fftw_malloc(sizeof(double)*windowSize);
  windowData= (double*) fftw_malloc(sizeof(double)*windowSize);
  plan_r2hc = fftw_plan_r2r_1d(windowSize, inputData, windowData, FFTW_R2HC, FFTW_PATIENT);
  // Opning files 
  inputFile = fopen("sinusD","rb");
  outputFile= fopen("windowingResult","wb+");
  if(inputFile==NULL ){
    printf("Couldn't open either the input or the output file \n");
    return -1;
  }

  while((i=fread(inputData,sizeof(double),windowSize,inputFile))==windowSize){   

    for( index1 =0; index1 < WINDOW_SIZE;index1++){
      inputData[index1]*=w[index1];
    //  printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     fftw_execute_r2r(plan_r2hc, inputData, windowData);
     for( index1 =0; index1 < windowSize;index1++){
          outputData[index1]+=windowData[index1];
     }
     if(overlaping!=0)
      fseek(inputFile,(-overlaping)*sizeof(double),SEEK_CUR);
  }
  if( i!=0){
     i = -i;
    fseek(inputFile ,i*sizeof(double),SEEK_END);
    fread(inputData,sizeof(double),-i,inputFile);
    for( index1 =0; index1 < -i;index1++){
      inputData[index1]*=(1.0 - cos(2.0 * M_PI * index1/(windowSize-1)) * 0.5);
    //  printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     fftw_execute_r2r(plan_r2hc, inputData, windowData);
     for( index1 =0; index1 < windowSize;index1++){
          outputData[index1]+=windowData[index1];
     }

  }


  powVal = outputData[0]*outputData[0];
  powVal /= (windowSize*windowSize)/2;
  index1 = 0;
  fprintf(outputFile,"%lf ",powVal);
  printf(" PSD \t %lf\n",powVal);
  for (index1 =1; index1<=windowSize;index1++){
            powVal = outputData[index1]*outputData[index1]+outputData[windowSize-index1]*outputData[windowSize- index1];
            powVal/=(windowSize*windowSize)/2;
          //  powVal = 20*log10(fabs(powVal));
            fprintf(outputFile,"%lf ",powVal);
            printf(" PsD %d \t %10.5lf\n",index1,powVal);
    }


  fftw_free(inputData);
  fftw_free(outputData);
  fftw_free(windowData);
  fclose(inputFile);
  fclose(outputFile);
}

为什么我得到的结果只有 0.0000？我问了很多关于窗口以及如何使用它的问题，但我真的不明白我在这里做错了什么，知道吗？

2.更新在 SleuthEye 的回答之后，结果看起来很不错，将结果与 MATLAB 使用的结果进行比较:

  [output,f] = pwelch(input,hann(8192));
    plot(output);

在导入 c 程序的结果后，PSD 相同但比例不同:

.

如您所见，比例不一样。

Answer 1

正如 chux 提到的，outputData 数组的元素没有被初始化。

此外，使用 Bartlett's method 获得的功率谱密度估计值应该平均频谱值的功率(而不是计算平均值的功率):

outputData[index1] += windowData[index1]*windowData[index1];

最后，如果频谱缩放对您的应用很重要(即，如果您需要的不仅仅是频率分量的相对强度)，那么您还应该将窗口效应考虑到归一化因子中，如 Numerical Recipes 中所述。 :

double Wss = 0.0;
for (i=0; i<WINDOW_SIZE; i++) {
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;

因此，将所有内容放在一起并考虑 Half-complex packing您使用的格式:

outputData = fftw_malloc(sizeof(double)*(windowSize/2 + 1));
for (index1 =0; index1 <= windowSize/2; index1++) {
  outputData[index1] = 0.0;
}
Wss = 0.0;
for (i=0; i<WINDOW_SIZE; i++) {
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;
...

count = 0;
while((i=fread(inputData,sizeof(double),windowSize,inputFile))==windowSize) {   
  for( index1 =0; index1 < WINDOW_SIZE;index1++){
    inputData[index1]*=w[index1];
  }
  fftw_execute_r2r(plan_r2hc, inputData, windowData);

  outputData[0] += windowData[0]*windowData[0];
  for( index1 =1; index1 < windowSize/2;index1++)
  {
    double re = windowData[index1];
    double im = windowData[windowSize-index1];
    outputData[index1] += re*re + im*im;
  }
  outputData[windowSize/2] += windowData[windowSize/2]*windowData[windowSize/2];
  count++;
}
...
if (halfSpectrum){
  norm = count*Wss/2;
  powVal = outputData[0]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
  for (index1 =1; index1<windowSize/2;index1++){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  powVal = outputData[windowSize/2]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
}
else{
  norm = count*Wss;
  for (index1 =0; index1<=windowSize/2;index1++){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  for (index1 =windowSize/2-1; index1>0;index1--){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
}

更新(解释与 Matlab pwelch 的差异):

根据 Octave 的 pwelch 帮助(应该与 Matlab 的输出匹配):

The spectral density is the mean of the periodograms, scaled so that area under the spectrum is the same as the mean square of the data.

换句话说，

而上面提供的比例因子适用于这样的定义，其中比例是这样的，即离散光谱值的总和与数据的均方相同(与this previous question的预期一致):

因此，定义上的差异引入了一个额外的 WINDOW_SIZE*sampling_rate 因子(请注意，如果不指定 pwelch 的第四个参数，您将使用默认的 sampling_rate 即 1Hz)。

因此，要使 C 版本的半谱输出匹配 pwelch，您需要:

norm = count*Wss/(2*WINDOW_SIZE*sampling_rate);
powVal = outputData[0]/(2*norm);
fprintf(outputFile,"%lf ",powVal);
for (index1 =1; index1<windowSize/2;index1++){
  powVal = outputData[index1]/norm;
  fprintf(outputFile,"%lf ",powVal);
}
powVal = outputData[windowSize/2]/(2*norm);
fprintf(outputFile,"%lf ",powVal);

否则，缩放 pwelch 输出以匹配 C 程序中使用的定义:

% For arbitrary sampling_rate:
%[output,f] = pwelch(input,hann(8192),[],[],sampling_rate)/(8192*sampling_rate);
% which simplifies to the following by setting sampling_rate = 1
[output,f] = pwelch(input,hann(8192))/8192;
plot(output);