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.
