实现开窗的正确方法

Question

我正在尝试在一个程序中实现窗口化，为此我编写了一个包含 2048 个样本的 sin 函数。我正在读取值并尝试使用“矩形”窗口计算 PSD，当我的窗口为 2048 宽时，结果是准确的。否则结果对我来说没有任何意义。

这是我正在使用的代码，

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


int main (){
  FILE* inputFile = NULL;
  FILE* outputFile= NULL;
  double* inputData=NULL; 
  double* outputData=NULL;
  double* windowData=NULL;
  unsigned int windowSize = 512;
  int overlaping =128;
  int index1 =0,index2=0, i=0;
  double powVal= 0.0;
  fftw_plan plan_r2hc;

// 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){   
      fftw_execute_r2r(plan_r2hc, inputData, windowData);
    for( index1 =0; index1 < windowSize;index1++){
      outputData[index1]+=windowData[index1];
      printf("index %d \t %lf\n",index1,inputData[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);


  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/2;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);
}

Answer 1

您需要用 window function 预乘信号.如果您正在计算多个 FFT，这可以预先计算。

例如，一个汉宁窗的计算方式如下:

#define WINDOW_SIZE 2048
int i;
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;
}

在计算傅立叶变换之前，按如下方式将输入数据乘以该窗口:

for (i=0; i<WINDOW_SIZE; i++) inputData[i] *= w[i];

解释

当您计算有限样本集的傅立叶变换时，您实际得到的是无限信号的频谱，您将永远重复这些样本。除非您采样的信号频率恰好是采样帧速率的倍数，否则您将在一个采样帧的末尾运行到下一个采样帧的开始处出现较大的不连续性。窗口函数使样本框边缘的样本变平，以消除这些不连续性。