Caculating le PSD en utilisant FFTW

Question

Je suis fichier audio que j'ai enregistrement en utilisant le `ALSA lib en utilisant les configurations suivantes:

Fs = 96000; // sample frequency 
channelNumber = 1 ; 
format =int16 ; 
length = 5sec; 

ce qui signifie que je reçois 480000 valeur 16bit. Maintenant, je veux calculer le PSD de l'ensemble de ce pour obtenir quelque chose comme:

ce que je suis en train de faire est TTO enregistrer le résultat en tant que groupe de valeur double en données supplémentaires, donc je peut les tracer de les évaluer (je ne sais pas si c'est exact):

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <fftw3.h> 

int main(){ 
    char fileName[] = "sound.raw"; 
    char magnFile[] = "data.txt"; 
    FILE* inp = NULL; 
    FILE* oup = NULL; 
    float* data = NULL; 
    fftwf_complex* out; 
    int index = 0; 
    fftwf_plan plan; 
    double var =0; 
    short wert = 0; 
    float r,i,magn; 
    int N = 512; 

    data =(float*)fftwf_malloc(sizeof(float)*N); 



    out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*N); 
    //Allocating the memory for the input data 
    plan = fftwf_plan_dft_r2c_1d(N,data,out, FFTW_MEASURE); 
    // opening the file for reading 
    inp = fopen(fileName,"r"); 
    oup = fopen(magnFile,"w+"); 

    if(inp== NULL){ 
     printf(" couldn't open the file \n "); 
     return -1; 
    } 
    if(oup==NULL){ 
     printf(" couldn't open the output file \n"); 
    } 
    while(!feof(inp)){ 

      if(index < N){ 
       fread(&wert,sizeof(short),1,inp); 
       //printf(" Wert %d \n",wert); 
       data[index] = (float)wert; 
       //printf(" Wert %lf \n",data[index]); 
       index = index +1; 
      } 
      else{ 

       index = 0; 
       fftwf_execute(plan); 
       //printf("New Plan \n"); 
       //printf(" Real \t imag \t Magn \t \n"); 
       for(index = 0 ; index<N; index++){ 
        r=out[index][0]; 
        i =out[index][1]; 
        magn = sqrt((r*r)+(i*i)); 
        printf("%.10lf \t %.10lf \t %.10lf \t \n",r,i,magn); 
        //fwrite(&magn,sizeof(float),1,oup); 
        //fwrite("\n",sizeof(char),1,oup); 
        fprintf(oup,"%.10lf\n ", magn); 
       } 
       index = 0 ; 
       fseek(inp,N,SEEK_CUR); 

      } 
    } 
    fftwf_destroy_plan(plan); 
    fftwf_free(data); 
    fftwf_free(out); 
    fclose(inp); 
    fclose(oup); 
    return 0 ; 
} 

le problème que j'ai est comment puis-je mettre en œuvre la fonction de remontage dans mon code? et je ne pense pas que ce résultat soit précis, puisque j'ai beaucoup de zéro dans les valeurs de magnitude? ?
Si quelqu'un a un exemple, je serai reconnaissant.

Answer 1

Voici un exemple simple d'appliquer un "Hanning" window à vos données avant la FFT:

for (int i = 0; i < N; ++i) 
{ 
    data[i] *= 0.5 * (1.0 + cos(2.0 * M_PI * (double)i/(double)(N - 1))); 
}