利用FFTW计算PSD

Question

我是使用“`ALSA”库使用以下设置记录的声音文件：

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

这意味着我得到了480000 16位的值。现在，我想计算出这组数据的PSD值，得到这样的结果：

​

​

我想要做的是将结果保存为额外数据中的一组双值，这样我就可以绘制它们的评估图(我不确定这是否正确)：

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

我遇到的问题是如何在代码中实现缠绕函数？我不认为这个结果是准确的，因为我得到了很多数量级的零值？？

如果有人有个榜样，我会感激的。

Answer 1

下面是将“汉宁”窗应用于您的数据的一个简单示例：

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