使用 FFTW 计算 PSD

Question

我是使用以下设置使用“ALSA”库录制的声音文件:

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

意思是我得到了 480000 个 16 位值。现在我想计算那组的 PSD 以获得类似的东西:

我想做的是将结果保存为额外数据中的一堆 double 值，这样我就可以绘制它们并对其进行评估(我不确定这是否正确):

   #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

这是一个应用 "Hanning" window 的简单示例FFT 之前的数据:

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