计算功率谱密度 [英] Calculating the Power spectral density

问题描述

我想获得一个真正的数据通过使用 fftw3 设定的PSD库
为了测试我写了如下图所示的小程序，产生了一个它遵循正弦函数信号

I am trying to get the PSD of a real data set by making use of fftw3 library To test I wrote a small program as shown below ,that generates the a signal which follows sinusoidal function

#include <stdio.h>
#include <math.h>
#define PI 3.14

int main (){
    double  value= 0.0;
    float frequency = 5;
    int i = 0 ; 
    double time = 0.0;
    FILE* outputFile = NULL;
    outputFile = fopen("sinvalues","wb+");
    if(outputFile==NULL){
        printf(" couldn't open the file \n");
        return -1;
    }

    for (i = 0; i<=5000;i++){

        value =  sin(2*PI*frequency*zeit);
        fwrite(&value,sizeof(double),1,outputFile);
        zeit += (1.0/frequency);
    }
    fclose(outputFile);
    return 0;

}

现在我在读以上程序的输出文件，并试图计算其PSD像如下图所示。

Now I'm reading the output file of above program and trying to calculate its PSD like as shown below

#include <stdio.h>
#include <fftw3.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14
int main (){
    FILE* inp = NULL;
    FILE* oup = NULL;
    double* value;// = 0.0;
    double* result;
    double spectr = 0.0 ;
    int windowsSize =512;
    double  power_spectrum = 0.0;
    fftw_plan plan;

    int index=0,i ,k;
    double multiplier =0.0;
    inp = fopen("1","rb");
    oup = fopen("psd","wb+");

    value=(double*)malloc(sizeof(double)*windowsSize);
    result = (double*)malloc(sizeof(double)*(windowsSize)); // what is the length that I have to choose here ? 
        plan =fftw_plan_r2r_1d(windowsSize,value,result,FFTW_R2HC,FFTW_ESTIMATE);

    while(!feof(inp)){

        index =fread(value,sizeof(double),windowsSize,inp);
            // zero padding 
        if( index != windowsSize){
            for(i=index;i<windowsSize;i++){
                    value[i] = 0.0;
                        }

        }


        // windowing  Hann 

        for (i=0; i<windowsSize; i++){
            multiplier = 0.5*(1-cos(2*PI*i/(windowsSize-1)));
            value[i] *= multiplier;
        }


        fftw_execute(plan);


        for(i = 0;i<(windowsSize/2 +1) ;i++){ //why only tell the half size of the window
            power_spectrum = result[i]*result[i] +result[windowsSize/2 +1 -i]*result[windowsSize/2 +1 -i];
            printf("%lf \t\t\t %d \n",power_spectrum,i);
            fprintf(oup," %lf \n ",power_spectrum);
        }

    }
    fclose(oup);
    fclose(inp);
    return 0;

}

荫不知道的，我这样做的方式是否正确，但低于是我所得到的结果：

Iam not sure about the correctness of the way I am doing this, but below are the results i have obtained:

任何一个可以帮助我在追查上述方法的误差

Can any one help me in tracing the errors of the above approach

在此先感谢
*更新
哈特穆特答案后I'vve编辑code，但仍然得到了相同的结果：

Thanks in advance *UPDATE after hartmut answer I'vve edited the code but still got the same result :

和输入数据如下：

更新
增加样本frequencyand 2048这里窗口的大小后，我得到了什么：

更新
使用ADD-ON这里经过怎样的结果看起来像使用窗口：

UPDATE after increasing the sample frequencyand a windows size of 2048 here is what I've got : UPDATE after using the ADD-ON here how the result looks like using the window :

推荐答案

您结合了错误的输出值功率谱线。有 windowsSize / 2 + 1 在结果开始实际值和 windowsSize / 2 - 在以相反的顺序结束1 虚值。这是因为，第一（0Hz时）的虚数分量和最后（奈奎斯特频率）频谱线是0

You combine the wrong output values to power spectrum lines. There are windowsSize / 2 + 1 real values at the beginning of result and windowsSize / 2 - 1 imaginary values at the end in reverse order. This is because the imaginary components of the first (0Hz) and last (Nyquist frequency) spectral lines are 0.

int spectrum_lines = windowsSize / 2 + 1;
power_spectrum = (double *)malloc( sizeof(double) * spectrum_lines );

power_spectrum[0] = result[0] * result[0];
for ( i = 1 ; i < windowsSize / 2 ; i++ )
    power_spectrum[i] = result[i]*result[i] + result[windowsSize-i]*result[windowsSize-i];
power_spectrum[i] = result[i] * result[i];

和有一个小的失误：你应该申请窗口功能只对输入信号而不是零填充部分

And there is a minor mistake: You should apply the window function only to the input signal and not to the zero-padding part.

ADD-ON：

您的测试程序产生正弦信号的样本5001，然后你阅读和分析这个信号的第512个样本。这样做的结果是，你只分析一个周期的一小部分。由于信号的硬截止它包含的能量几乎未predictable能级范围广泛，因为你甚至不会使用PI，但只有3.41而不是precise够做任何predictable计算。

Your test program generates 5001 samples of a sinusoid signal and then you read and analyse the first 512 samples of this signal. The result of this is that you analyse only a fraction of a period. Due to the hard cut-off of the signal it contains a wide spectrum of energy with almost unpredictable energy levels, because you not even use PI but only 3.41 which is not precise enough to do any predictable calculation.

您需要保证周期的整数数目恰好装配到512个样本的分析窗口。因此，你应该在你的测试信号创建程序更改为具有完全相同 numberOfPeriods 在测试信号周期（如 numberOfPeriods = 1 意味着sinoid的一个周期的周期正好512个样本，2 => 256，3 =>3分之512，4 => 128，...）的。这样一来，您就可以在特定的光谱线来产生能量。请记住， windowSize 必须在这两个程序相同的值，因为不同的尺寸做这方面的努力也没用。

You need to guarantee that an integer number of periods is exactly fitting into your analysis window of 512 samples. Therefore, you should change this in your test signal creation program to have exactly numberOfPeriods periods in your test signal (e.g. numberOfPeriods=1 means that one period of the sinoid has a period of exactly 512 samples, 2 => 256, 3 => 512/3, 4 => 128, ...). This way, you are able to generate energy at a specific spectral line. Keep in mind that windowSize must have the same value in both programs because different sizes make this effort useless.

#define PI 3.141592653589793 // This has to be absolutely exact!

int windowSize = 512;        // Total number of created samples in the test signal
int numberOfPeriods = 64;    // Total number of sinoid periods in the test signal
for ( n = 0 ; n < windowSize ; ++n ) {
    value = sin( (2 * PI * numberOfPeriods * n) / windowSize );
    fwrite( &value, sizeof(double), 1, outputFile );
}

这篇关于计算功率谱密度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！