# Calculating the Power spectral density [closed]

I've already asked this question on SO, but didn't get an answer, I tqake a shot here. 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;
float  length = 2;
int index = 4;
double 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)){

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

and the input data look like :

• In your code listing the sampling frequency and the frequency of your tone are not defined very well. zeit is not initialized. Also zeit = n/F, where n is an integer. Then $\sin(2\pi*F*zeit)= \sin(2\pi*F*n/F)=\sin(2\pi*n)=0$, so the output will be a constant. Commented Jan 7, 2015 at 20:38
• your sine production is essentailly broken, as @David says. $\sin(2\pi f \frac nf)\equiv \sin(2\pi n)\equiv 0 \quad \forall n\in\mathbb N$. the fact you're seeing something like a sine at all is only due to the limited machine resolution! Commented Jul 31, 2016 at 14:56
• I'm voting to close this question as off-topic because the data analyzed here is produced in a broken way and answering the question hence doesn't solve any relevant problem. Commented Jul 31, 2016 at 14:57