使用 FFTW 计算 PSD

Question

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

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

意思是我得到 480000 16bit 值。现在我想计算该集合的 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

以下是在 FFT 之前将"Hanning" window 应用于您的数据的简单示例：

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