Estou tentando obter o PSD de um conjunto de dados real usandofftw3 library Para testar, escrevi um pequeno programa como mostrado abaixo, que gera o sinal a que segue a função sinusoidal

#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;

}

Agora estou lendo o arquivo de saída do programa acima e tentando calcular seu PSD como mostrado abaixo

#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;

}

Não tenho certeza sobre a correção da maneira como estou fazendo isso, mas abaixo estão os resultados que obtive:

Alguém pode me ajudar a rastrear os erros da abordagem acima

desde já, obrigado*ATUALIZAR depois da resposta do hartmut, editei o código, mas ainda assim obtive o mesmo resultado:

e os dados de entrada se parecem com:

ATUALIZAR depois de aumentar a frequência da amostra e um tamanho de janela de 2048, aqui está o que eu tenho: ATUALIZAR depois de usar o ADD-ON aqui, como é o resultado usando a janela: