PSD using FFTW Halfcomplex transformation

Question

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;
}

the value that I get are :

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

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

thanks in advance !

UPDATE

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]);
        }

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 !

score 2 · Accepted Answer · answered Aug 08 '14 at 02:50

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:
enter image description here

where P(k) is related to the FFT outputs X(k) through:
enter image description here
for some scaling factor S.

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

Applying Parseval's theorem:
enter image description here
to the definition yields:

Or S=1/N².

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.

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);

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).

thanks for your answer this and the other one, need a bit time to understand it — Engine, Aug 08 '14 at 10:51

PSD using FFTW Halfcomplex transformation

1 Answers1

Linked