Calculating the Power spectral density

Question

I am trying to get the PSD of a real data set by making use of fftw3 library To test I wrote a small program as shown below ,that generates the a signal which follows sinusoidal function

#include <stdio.h>
#include <math.h>
#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<=5000;i++){

        value =  sin(2*PI*frequency*zeit);
        fwrite(&value,sizeof(double),1,outputFile);
        zeit += (1.0/frequency);
    }
    fclose(outputFile);
    return 0;

}

Now I'm reading the output file of above program and trying to calculate its PSD like as shown below

#include <stdio.h>
#include <fftw3.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#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<windowsSize;i++){
                    value[i] = 0.0;
                        }

        }


        // windowing  Hann 

        for (i=0; i<windowsSize; i++){
            multiplier = 0.5*(1-cos(2*PI*i/(windowsSize-1)));
            value[i] *= multiplier;
        }


        fftw_execute(plan);


        for(i = 0;i<(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;

}

Iam not sure about the correctness of the way I am doing this, but below are the results i have obtained:

Can any one help me in tracing the errors of the above approach

Thanks in advance *UPDATE after hartmut answer I'vve edited the code but still got the same result :

and the input data look like :

UPDATE after increasing the sample frequencyand a windows size of 2048 here is what I've got : UPDATE after using the ADD-ON here how the result looks like using the window :

Answer 1

You combine the wrong output values to power spectrum lines. There are windowsSize / 2 + 1 real values at the beginning of result and windowsSize / 2 - 1 imaginary values at the end in reverse order. This is because the imaginary components of the first (0Hz) and last (Nyquist frequency) spectral lines are 0.

int spectrum_lines = windowsSize / 2 + 1;
power_spectrum = (double *)malloc( sizeof(double) * spectrum_lines );

power_spectrum[0] = result[0] * result[0];
for ( i = 1 ; i < windowsSize / 2 ; i++ )
    power_spectrum[i] = result[i]*result[i] + result[windowsSize-i]*result[windowsSize-i];
power_spectrum[i] = result[i] * result[i];

And there is a minor mistake: You should apply the window function only to the input signal and not to the zero-padding part.

ADD-ON:

Your test program generates 5001 samples of a sinusoid signal and then you read and analyse the first 512 samples of this signal. The result of this is that you analyse only a fraction of a period. Due to the hard cut-off of the signal it contains a wide spectrum of energy with almost unpredictable energy levels, because you not even use PI but only 3.41 which is not precise enough to do any predictable calculation.

You need to guarantee that an integer number of periods is exactly fitting into your analysis window of 512 samples. Therefore, you should change this in your test signal creation program to have exactly numberOfPeriods periods in your test signal (e.g. numberOfPeriods=1 means that one period of the sinoid has a period of exactly 512 samples, 2 => 256, 3 => 512/3, 4 => 128, ...). This way, you are able to generate energy at a specific spectral line. Keep in mind that windowSize must have the same value in both programs because different sizes make this effort useless.

#define PI 3.141592653589793 // This has to be absolutely exact!

int windowSize = 512;        // Total number of created samples in the test signal
int numberOfPeriods = 64;    // Total number of sinoid periods in the test signal
for ( n = 0 ; n < windowSize ; ++n ) {
    value = sin( (2 * PI * numberOfPeriods * n) / windowSize );
    fwrite( &value, sizeof(double), 1, outputFile );
}