C语言 实现窗口的正确方法

Question

我试图在程序中实现窗口，为此我编写了一个带有2048个样本的sin函数。我正在读取值并尝试使用"rect"窗口计算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 %lfn",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 %lfn",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.5lfn",index1,powVal);
    }

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

Answer 1

需要用窗函数对信号进行预乘。如果您正在计算多个fft，则可以预先计算此值。

例如，Hanning窗口的计算方法如下:

#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];

当你计算一组有限样本的傅里叶变换时，你实际上得到的是无限信号的频谱你将通过永远重复这些样本得到这些频谱。除非你采样的信号的频率是采样帧率的精确倍数，否则你会得到很大的不连续，一个采样帧的结束会进入下一个采样帧的开始。窗函数使样本帧边缘的样本变平，以消除这些不连续。