title: Calculating the Power spectral density
body: <p>I am trying to get the PSD of a real data set by making use of  <code>fftw3</code> library 
  To test I wrote a small program as shown below ,that generates the  a signal which follows sinusoidal function </p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;math.h&gt;
  #define PI 3.14
  
  int main (){
      double  value= 0.0;
      float frequency = 5;
      int i = 0 ; 
      double time = 0.0;
      FILE* outputFile = NULL;
      outputFile = fopen("sinvalues","wb+");
      if(outputFile==NULL){
          printf(" couldn't open the file \n");
          return -1;
      }
  
      for (i = 0; i&lt;=5000;i++){
  
          value =  sin(2*PI*frequency*zeit);
          fwrite(&amp;value,sizeof(double),1,outputFile);
          zeit += (1.0/frequency);
      }
      fclose(outputFile);
      return 0;
  
  }
  </code></pre>
  
  <p>Now I'm reading the output file of above program and trying to calculate its PSD like as shown below </p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;fftw3.h&gt;
  #include &lt;complex.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;math.h&gt;
  #define PI 3.14
  int main (){
      FILE* inp = NULL;
      FILE* oup = NULL;
      double* value;// = 0.0;
      double* result;
      double spectr = 0.0 ;
      int windowsSize =512;
      double  power_spectrum = 0.0;
      fftw_plan plan;
  
      int index=0,i ,k;
      double multiplier =0.0;
      inp = fopen("1","rb");
      oup = fopen("psd","wb+");
  
      value=(double*)malloc(sizeof(double)*windowsSize);
      result = (double*)malloc(sizeof(double)*(windowsSize)); // what is the length that I have to choose here ? 
          plan =fftw_plan_r2r_1d(windowsSize,value,result,FFTW_R2HC,FFTW_ESTIMATE);
  
      while(!feof(inp)){
  
          index =fread(value,sizeof(double),windowsSize,inp);
              // zero padding 
          if( index != windowsSize){
              for(i=index;i&lt;windowsSize;i++){
                      value[i] = 0.0;
                          }
  
          }
  
  
          // windowing  Hann 
  
          for (i=0; i&lt;windowsSize; i++){
              multiplier = 0.5*(1-cos(2*PI*i/(windowsSize-1)));
              value[i] *= multiplier;
          }
  
  
          fftw_execute(plan);
  
  
          for(i = 0;i&lt;(windowsSize/2 +1) ;i++){ //why only tell the half size of the window
              power_spectrum = result[i]*result[i] +result[windowsSize/2 +1 -i]*result[windowsSize/2 +1 -i];
              printf("%lf \t\t\t %d \n",power_spectrum,i);
              fprintf(oup," %lf \n ",power_spectrum);
          }
  
      }
      fclose(oup);
      fclose(inp);
      return 0;
  
  }
  </code></pre>
  
  <p>Iam not sure about the correctness of the way I am doing this, but below are the results i have obtained:</p>
  
  <p><img src="https://i.stack.imgur.com/5VBKI.png" alt="PSD"></p>
  
  <p>Can any one help me in tracing the errors of the above approach </p>
  
  <p>Thanks in advance 
  <strong>*UPDATE</strong>
  after hartmut answer  I'vve edited the code but still got the same result :
  <img src="https://i.stack.imgur.com/IItkx.png" alt="PSD2"></p>
  
  <p>and the input data look like : </p>
  
  <p><img src="https://i.stack.imgur.com/penpP.png" alt="sinus"></p>
  
  <p><strong>UPDATE</strong>
  after  increasing the sample frequencyand a windows size of 2048 here is what I've got :
  <img src="https://i.stack.imgur.com/bR18q.png" alt="enter image description here">
  <strong>UPDATE</strong>
  after  using the ADD-ON here how the result looks like using the window :
  <img src="https://i.stack.imgur.com/v0bpO.png" alt="enter image description here"></p>