Need help making this Program calculating PSD

Question

I'm trying to make a program that calculate a PSD of a time serie(16384 sample), say a sinus here is the code :

// generating sin samples 

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

int main(){
    FILE *inp =NULL,*inp2;
    double value = 0.0; 
    float frequency = 1; // signal frequency 
    double timeSec = 1 ; // time in sec 
    unsigned int  numberOFSamples = 2048*8;
    double steps = timeSec / numberOFSamples; 
    double timer =0.0;
    float  dcValue =0.0;
    double index = 0;
     inp = fopen("sinus","wb+");
     inp2= fopen("sinusD","wb+");
    for( timer=0.0 ; timer<timeSec;timer+=steps){
        value= sin(2*M_PI*frequency*timer) +dcValue;
        fprintf(inp,"%lf ",value);
        fwrite(&value,sizeof(double),1,inp2);

     }

     fclose(inp);
     fclose(inp2);
     return 0;

}

the generating of the sinus values is correct I've check it using Matlab now to the PSD should be 1024 large here is the code :

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

#define WINDOW_SIZE 1024

int main (){
  FILE* inputFile = NULL;
  FILE* outputFile= NULL;
  double* inputData=NULL; 
  double* outputData=NULL;
  double* windowData=NULL;
  unsigned int windowSize = 1024;
  int overlaping =0;
  int index1 =0,index2=0, i=0;
  double powVal= 0.0;
  fftw_plan plan_r2hc;
  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;  // hann window
  }

// 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){   

    for( index1 =0; index1 < WINDOW_SIZE;index1++){
      inputData[index1]*=w[index1];
    //  printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     fftw_execute_r2r(plan_r2hc, inputData, windowData);
     for( index1 =0; index1 < windowSize;index1++){
          outputData[index1]+=windowData[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);
    for( index1 =0; index1 < -i;index1++){
      inputData[index1]*=(1.0 - cos(2.0 * M_PI * index1/(windowSize-1)) * 0.5);
    //  printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     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;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);
}

the result that I get is only 0.0000 why ? I'Ve asked mutiple question about window and how to use it, but I really don't get what I'm doing wrong here , any idea ?

2.UPDATE After SleuthEye's answer the result looks pretty good, comparing the result and the the result of MATLAB using :

  [output,f] = pwelch(input,hann(8192));
    plot(output);

and after import the result of c program the PSD is the same but the scale is different :

.

As you can see the scale isn't the same.

Answer 1

As chux mentionned, elements of the outputData array are not initialized.

Furthermore, power-spectrum density estimates obtained using Bartlett's method should average the power of the spectrum values (rather than computing the power of the average):

outputData[index1] += windowData[index1]*windowData[index1];

Finally, if spectrum scaling is important for your application (ie. if you need more than the relative strength of the frequency components), then you should also take into account windowing effect into your normalization factor, as described in Numerical Recipes:

double Wss = 0.0;
for (i=0; i<WINDOW_SIZE; i++) {
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;

So, putting it all together and taking into account the Half-complex packing format which you are using:

outputData = fftw_malloc(sizeof(double)*(windowSize/2 + 1));
for (index1 =0; index1 <= windowSize/2; index1++) {
  outputData[index1] = 0.0;
}
Wss = 0.0;
for (i=0; i<WINDOW_SIZE; i++) {
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;
...

count = 0;
while((i=fread(inputData,sizeof(double),windowSize,inputFile))==windowSize) {   
  for( index1 =0; index1 < WINDOW_SIZE;index1++){
    inputData[index1]*=w[index1];
  }
  fftw_execute_r2r(plan_r2hc, inputData, windowData);

  outputData[0] += windowData[0]*windowData[0];
  for( index1 =1; index1 < windowSize/2;index1++)
  {
    double re = windowData[index1];
    double im = windowData[windowSize-index1];
    outputData[index1] += re*re + im*im;
  }
  outputData[windowSize/2] += windowData[windowSize/2]*windowData[windowSize/2];
  count++;
}
...
if (halfSpectrum){
  norm = count*Wss/2;
  powVal = outputData[0]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
  for (index1 =1; index1<windowSize/2;index1++){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  powVal = outputData[windowSize/2]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
}
else{
  norm = count*Wss;
  for (index1 =0; index1<=windowSize/2;index1++){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  for (index1 =windowSize/2-1; index1>0;index1--){
    powVal = outputData[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
}

Update (explanation of discrepancy with Matlab pwelch):

According to pwelch help from Octave (which should match Matlab's output):

The spectral density is the mean of the periodograms, scaled so that area under the spectrum is the same as the mean square of the data.

In other words,

Whereas the scaling factor provided above applies for a definition where the scaling is such that the sum of discrete spectrum values the same as the mean square of the data (consistently with the expectation of this previous question):

Thus, the difference in definition introduces an extra WINDOW_SIZE*sampling_rate factor (note that without specifying the fourth argument to pwelch, you are using the default sampling_rate which is 1Hz).

So, for the half spectrum output of the C version to match pwelch you would need:

norm = count*Wss/(2*WINDOW_SIZE*sampling_rate);
powVal = outputData[0]/(2*norm);
fprintf(outputFile,"%lf ",powVal);
for (index1 =1; index1<windowSize/2;index1++){
  powVal = outputData[index1]/norm;
  fprintf(outputFile,"%lf ",powVal);
}
powVal = outputData[windowSize/2]/(2*norm);
fprintf(outputFile,"%lf ",powVal);

Or else, to scale pwelch output to match the definition used in the C program:

% For arbitrary sampling_rate:
%[output,f] = pwelch(input,hann(8192),[],[],sampling_rate)/(8192*sampling_rate);
% which simplifies to the following by setting sampling_rate = 1
[output,f] = pwelch(input,hann(8192))/8192;
plot(output);

Answer 2

Elements of outputData[] are not initialized.

outputData = (double*) fftw_malloc(sizeof(double)*windowSize);
...
  outputData[index1] += ...

Suggest:

outputData = fftw_malloc(sizeof(double)*windowSize);
for (index1 =0; index1 < windowSize; index1++) {
  outputData[index1] = 0.0;
}
...
  outputData[index1] += ...