PSD使用FFTW Halfcomplex转型 [英] PSD using FFTW Halfcomplex transformation

问题描述

我问一个<一个href=\"http://stackoverflow.com/questions/24696122/calculating-the-power-spectral-density/24739037?noredirect=1#comment38469909_24739037\">similar问题，这是回答了，但是当我尝试做我的方式，我得到奇怪的值。
我希望得到一个正弦波的PSD使用半复杂的转换，如：

I've asked a similar question,which was answered but when I try to do it my way I get "strange" values. I want to get the PSD of a sin wave use the half complex transformation like :

    #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值是零，峰值的位置是正确的，但他们的价值没有任何想法，为什么？

otherwise the PSD value is zero, the position of the peak is correct but their value isn't any idea why ?

在此先感谢！

更新

我解决这个问题，但我不知道为什么它的工作原理，所以我不会把它作为一个答案：
这里是我的code已经改变了：

I solve it but I don't get why it works , so I won't put it as an answer : here is what I've changed in the code :

.......................................
      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的广场和2分的解释吗？再次感谢您的帮助！

the result is accurate , I hope getting any explanation about why I should divide on the square of the windowsSize and the on 2 ? thanks again for your help !

推荐答案

如数字食谱描述，正常化功率谱密度（PSD）的可能因PSD你的确切定义有所不同。
一种可能的定义，这使得频谱估计值的总和对应于时域函数的均值的平方幅度：搜索

As described in Numerical Recipes, the normalization of the power-spectrum density (PSD) may vary depending on your exact definition of PSD. One possible definition which makes the sum of spectrum estimates correspond to the mean squared amplitude of the time domain function:

其中，P（k）的通过有关的FFT输出X（K）：结果，
结果
对于一些缩放系数取值。

where P(k) is related to the FFT outputs X(k) through:

for some scaling factor S.

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

This seems to be the definition you are using based on your expected results.

应用<一href=\"http://en.wikipedia.org/wiki/Discrete_Fourier_transform#The_Plancherel_theorem_and_Parseval.27s_theorem\"相对=nofollow> Parseval定理：结果
结果
的定义，得出：搜索
结果
或S = 1 / N ² 。

这当然假设你与整个频谱工作。
在另一方面 FFTW的半复杂的格式正如其名称所暗示的，给你只有一半的光谱
（在另一半对称的实值输入）。
通过在该对称频谱值的等值平方量值加入得到的频率分量的总功率。
这是你在​​哪里得到的2从因素。
请注意，有没有相应的实值指数0和的NumberOfSamples / 2 对称的输出，因此你应该2在这些情况下不繁殖。

This of course assumes that you work with the entire spectrum. FFTW's Half-complex format on the other hand, as the name implies, gives you only half the spectrum (the other half being symmetric for real-valued inputs). You get the total power of a frequency component by adding in the equal-valued squared-magnitude of that symmetric spectrum value. This is where you get the factor of 2 from. Note that there are no corresponding symmetric outputs for real-valued index 0 and numberOfSamples/2, so you should not multiply by 2 in these cases.

那么，你的PSD应计算为：

So, your PSD should be computed as:

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）

Note: the mean squared amplitude of your DC component should be 0.2*0.2 = 0.04 (not 0.2 which you indicated as your expected value for index 0).

这篇关于PSD使用FFTW Halfcomplex转型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！