Search code examples
csignal-processingfftw

Caculating the PSD using FFTW

I'm sound file that I've record using the `ALSA' lib using following setups :

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

meaning that I get 480000 16bit value. Now I want to calculate the PSD of the set of that to get something like :

PSD

what I'm trying to do is tto save the result as a bunch of double value in a extra data so I can plot them of evaluating them ( I'm not sure if that's correct) : 

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

the problem that I have is how can I implement the winding function in my code ? and I don't think that result is accurate, since I'get a lot of zero in magnitude values ? ?
if somebody has an example I'll be thankful .

Solution

  • Here is a simple example of applying a "Hanning" window to your data prior to the FFT:

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