使用 FFTW 半复变换的 PSD

Question

我问了一个similar 的问题，这个问题得到了回答，但是当我尝试按照自己的方式去做时，我得到了“奇怪”的值。 我想使用半复数变换来获得正弦波的 PSD：

    #include <stdio.h>
#include <fftw3.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793
int main (){

        double* inputValues;
        double* outputValues;
        double  realVal;
        double  imagVal;
        double  powVal=0.0;
        double  absVal;
        double timer;
        fftw_plan plan;
        double timeIntervall= 1.0; // 1sec 
        int numberOfSamples  =512;
        double timeSteps = timeIntervall/numberOfSamples;
        float frequency=10.0;
        float dcValue = 0.2;
        float value=0.0;
        int index=0;
        // allocating the memory for the fftw arrays 
        inputValues = (double*) fftw_malloc(sizeof(double)* numberOfSamples);
        outputValues = (double *) fftw_malloc(sizeof(double)*(numberOfSamples/*2*/));
        plan = fftw_plan_r2r_1d(numberOfSamples,inputValues,outputValues,FFTW_R2HC,FFTW_ESTIMATE);


    for (timer = 0; timer<=timeIntervall; timer += timeSteps){
        value =  sin(2*M_PI*frequency*timer) +dcValue;
        inputValues[index++] = value;

    }


        fftw_execute(plan);

        for (index=0;index<=numberOfSamples/*2*/;index++){
            powVal = outputValues[index]*outputValues[index]+outputValues[numberOfSamples-index]*outputValues[numberOfSamples- index];
            if(index==0)
                powVal/=2;
            powVal/=numberOfSamples;
            fprintf(stdout," index %d \t PSD value %lf \n",index,powVal);
        }

    return 0;
}

我得到的价值是：

index 0          PSD value 12.24  // expecting 0.2
................
.....................
index 10         PSD value 129.99999  // expecting 0.5
........
.......
index 502       PSD value 127.9999  // expecting 0.5
......................
......................
index 512       PSD value 12.24   // expecting 0.2

否则 PSD 值为零，峰值的位置是正确的，但它们的值不知道为什么？

提前致谢！

更新

我解决了它，但我不明白它为什么起作用，所以我不会把它作为答案： 这是我在代码中所做的更改：

.......................................
      fftw_execute_r2r(plan_r2hc, in, out);
  powVal = outputValues[0]*outputValues[0];
  powVal /= (numberOfSamples*numberOfSamples)/2;  ///WHY ??????
  index = 0;
fprintf(stdout," index %d \t PSD value %lf \t \t %lf \n",index,powVal,outputValues[index]);
  for (index =1; index<numberOfSamples/2;index++){
  powVal = outputValues[index]*outputValues[index]+outputValues[numberOfSamples-index]*outputValues[numberOfSamples- index];

            powVal/=(numberOfSamples*numberOfSamples)/2;  //WHY?????
            fprintf(stdout," index %d \t PSD value %lf \t \t %lf \n",index,powVal,outputValues[index]);
        }

结果是准确的，我希望得到任何解释，为什么我应该在 windowsSize 的平方和 on 2 上划分？再次感谢您的帮助 ！

Answer 1

如Numerical Recipes 中所述，功率谱密度 (PSD) 的归一化可能会因您对 PSD 的确切定义而异。 使频谱估计之和对应于时域函数的均方幅度的一种可能定义：

其中 P(k) 通过以下方式与 FFT 输出 X(k) 相关：

对于一些比例因子S。

这似乎是您根据预期结果使用的定义。

申请Parseval's theorem:

到定义产生：

或 S=1/N²。

这当然假设您使用整个频谱。 另一方面，FFTW's Half-complex format，顾名思义，只给你一半的频谱 （另一半对于实值输入是对称的）。 您可以通过添加该对称频谱值的等值平方幅度来获得频率分量的总功率。 这就是你得到因子 2 的地方。 请注意，实值索引 0 和 numberOfSamples/2 没有对应的对称输出，因此在这些情况下不应乘以 2。

因此，您的 PSD 应计算为：

powVal = outputValues[0]*outputValues[0];
powVal /= (numberOfSamples*numberOfSamples);
for (index =1; index<numberOfSamples/2;index++){
  powVal = outputValues[index]*outputValues[index]+outputValues[numberOfSamples-index]*outputValues[numberOfSamples-index];
  powVal /= (numberOfSamples*numberOfSamples)/2;
}
powVal = outputValues[numberOfSamples/2]*outputValues[numberOfSamples/2];
powVal /= (numberOfSamples*numberOfSamples);

注意：直流分量的均方幅度应为 0.2*0.2 = 0.04（而不是您指出的指数 0 的预期值的 0.2）。