计算功率谱密度

Question

我正在尝试通过使用 fftw3 库来获取真实数据集的 PSD 为了测试我写了一个如下所示的小程序，它生成一个遵循正弦函数的信号

#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，如下所示

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

}

我不确定我这样做的方式是否正确，但以下是我得到的结果：

任何人都可以帮助我追踪上述方法的错误

提前致谢 *更新 在 hartmut 回答之后，我已经编辑了代码，但仍然得到了相同的结果：

输入数据如下：

更新 在增加采样频率和 2048 的窗口大小之后，我得到了： 更新 在此处使用 ADD-ON 后，使用窗口的结果如何：

Answer 1

您将错误的输出值组合到功率谱线。 result 开头有windowsSize / 2 + 1 实数值，末尾有windowsSize / 2 - 1 虚数值，倒序排列。这是因为第一条（0Hz）和最后一条（奈奎斯特频率）谱线的虚部都是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];

还有一个小错误：您应该只将窗口函数应用于输入信号而不是零填充部分。

附加功能：

您的测试程序生成 5001 个正弦信号样本，然后您读取并分析该信号的前 512 个样本。这样做的结果是您只分析了一段时期的一小部分。由于信号的硬截止，它包含具有几乎不可预测的能量水平的广泛能量，因为您甚至不使用 PI，而仅使用 3.41，这不足以进行任何可预测的计算。

您需要保证整数个周期正好适合您的 512 个样本的分析窗口。因此，您应该在您的测试信号创建程序中将其更改为在您的测试信号中有准确的numberOfPeriods 周期（例如numberOfPeriods=1 表示正弦波的一个周期具有正好 512 个样本的周期，2 => 256, 3 = > 512/3, 4 => 128, ...)。这样，您就可以在特定的谱线上产生能量。请记住，windowSize 在两个程序中必须具有相同的值，因为不同的大小会使这项工作毫无用处。

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

Answer 2

对您预期的输出函数的一些说明。

• 您的输入是一个具有纯实数值的函数。 DFT 的结果具有复数值。 因此，您必须将变量 out 声明为不是 double 而是 fftw_complex *out。

• 通常，dft 输入值的数量与输出值的数量相同。 但是，dft 的输出频谱包含正的复振幅 频率以及负频率。

• 在纯实数输入的特殊情况下，正频率的幅度为 负频率幅度的共轭复数值。 为此，仅计算正频谱的频率， 这意味着复数输出值的数量是一半 实际输入值的数量。

• 如果您的输入是简单的正弦波，则频谱仅包含单个频率分量。 这适用于 10、100、1000 甚至更多的输入样本。 所有其他值为零。因此，处理大量输入值没有任何意义。

• 如果输入数据集包含单个周期，则复数输出值为 包含在 out[1] 中。

• 如果如果输入数据集包含 M 个完整的周期，在您的情况下为 5， 所以结果存储在 out[5]

• 我对您的代码做了一些修改。让一些事实更清楚。

---

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include "fftw3.h"

int performDFT(int nbrOfInputSamples, char *fileName)
{
    int nbrOfOutputSamples;
    double *in;
    fftw_complex *out;
    fftw_plan p;

    // In the case of pure real input data,
    // the output values of the positive frequencies and the negative frequencies
    // are conjugated complex values.
    // This means, that there no need for calculating both.
    // If you have the complex values for the positive frequencies,
    // you can calculate the values of the negative frequencies just by
    // changing the sign of the value's imaginary part
    // So the number of complex output values ( amplitudes of frequency components)
    // are the half of the number of the real input values ( amplitutes in time domain):
    nbrOfOutputSamples = ceil(nbrOfInputSamples/2.0);

    // Create a plan for a 1D DFT with real input and complex output
    in = (double*) fftw_malloc(sizeof(double) * nbrOfInputSamples);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nbrOfOutputSamples);
    p = fftw_plan_dft_r2c_1d(nbrOfInputSamples, in, out, FFTW_ESTIMATE);

    // Read data from input file to input array
    FILE* inputFile = NULL;
    inputFile = fopen(fileName,"r");
    if(inputFile==NULL){
        fprintf(stdout,"couldn't open the file %s\n", fileName);
        return -1;
    }
    double value;
    int     idx = 0;
    while(!feof(inputFile)){
        fscanf(inputFile, "%lf", &value);
        in[idx++] = value;       
    }
    fclose(inputFile);

    // Perform the dft
    fftw_execute(p);

    // Print output results
    char outputFileName[] = "dftvalues.txt";
    FILE* outputFile = NULL;
   outputFile = fopen(outputFileName,"w+");
    if(outputFile==NULL){
        fprintf(stdout,"couldn't open the file %s\n", outputFileName);
        return -1;
    }
    double  realVal;
    double  imagVal;
    double  powVal;
    double  absVal;

    fprintf(stdout, "     Frequency  Real       Imag        Abs       Power\n");
    for (idx=0; idx<nbrOfOutputSamples; idx++) {
        realVal = out[idx][0]/nbrOfInputSamples; // Ideed nbrOfInputSamples is correct!
        imagVal = out[idx][1]/nbrOfInputSamples; // Ideed nbrOfInputSamples is correct!
        powVal  = 2*(realVal*realVal + imagVal*imagVal);
        absVal  = sqrt(powVal/2);
        if (idx == 0) {
            powVal /=2;
        }
        fprintf(outputFile, "%10i %10.4lf %10.4lf %10.4lf %10.4lf\n", idx, realVal, imagVal, absVal, powVal);
        fprintf(stdout,     "%10i %10.4lf %10.4lf %10.4lf %10.4lf\n", idx, realVal, imagVal, absVal, powVal);
        // The total signal power of a frequency is the sum of the power of the posive and the negative frequency line.
        // Because only the positive spectrum is calculated, the power is multiplied by two.
        // However, there is only one single line in the prectrum for DC.
        // This means, the DC value must not be doubled.

    }
    fclose(outputFile);


    // Clean up
    fftw_destroy_plan(p);
    fftw_free(in); fftw_free(out);

    return 0;

}

int main(int argc, const char * argv[]) {
    // Set basic parameters
    float   timeIntervall = 1.0;     // in seconds
    int     nbrOfSamples = 50;      //  number of Samples per time intervall, so the unit is S/s
    double  timeStep = timeIntervall/nbrOfSamples; // in seconds
    float   frequency = 5;          // frequency in Hz

    // The period time of the signal is 1/5Hz = 0.2s
    // The number of samples per period is: nbrOfSamples/frequency = (50S/s)/5Hz = 10S
    // The number of periods per time intervall is: frequency*timeIntervall = 5Hz*1.0s = (5/s)*1.0s = 5



    // Open file for writing signal values
    char fileName[] = "sinvalues.txt";
    FILE* outputFile = NULL;
    outputFile = fopen(fileName,"w+");
    if(outputFile==NULL){
        fprintf(stdout,"couldn't open the file %s\n", fileName);
        return -1;
    }

    // Calculate signal values and write them to file
    double  time;
    double  value;
    double  dcValue = 0.2;
    int idx = 0;
    fprintf(stdout, "     SampleNbr   Signal value\n");
    for (time = 0; time<=timeIntervall; time += timeStep){
        value =  sin(2*M_PI*frequency*time) + dcValue;
        fprintf(outputFile, "%lf\n",value);
        fprintf(stdout, "%10i %15.5f\n",idx++, value);
    }
    fclose(outputFile);

    performDFT(nbrOfSamples, fileName);
    return 0;

}

• 如果 dft 的输入是纯实数，则输出无论如何都是复数。 所以你必须使用计划 r2c (RealToComplex)。
• 如果信号是sin(2*pi*f*t)，从t=0开始，频谱包含单条频率线 在 f，这是纯虚数。
• 如果符号有相位偏移，例如 sin(2*pi*f*t+phi)，则单行的值是复数。
• 如果您的采样频率为 fs，则输出频谱的范围为 -fs/2 ... +fs/2。
• 正负频率的实部相同。
• 正负频率的虚部符号相反。 这称为共轭复合物。
• 如果您有正谱的复数值，您可以计算 通过改变虚部的符号来获得负谱。 由于这个原因，不需要同时计算正谱和负谱。
• 一个边带保存所有信息。 因此，计划 r2c 中输出样本的数量是数量的一半+1 输入样本。

• 要获得频率的功率，您还必须考虑正频率 作为负频率。然而，计划 r2c 只提供正确的正半部分 的频谱。所以你必须将正极的功率加倍才能得到总功率。

  顺便说一句，fftw3 包的文档很好地描述了计划的用法。 您应该花时间阅读手册。

Answer 3

我不确定你的问题是什么。根据提供的信息，您的结果似乎是合理的。

您必须知道，PSD 是自相关函数的傅里叶变换。使用正弦波输入，您的 AC 函数将是周期性的，因此 PSD 会有音调，就像您绘制的那样。

我的“答案”确实是一些关于调试的想法。如果我们可以发布方程式，所有相关人员都会更容易。您可能知道这些天在 SE 上有一个信号处理部分。

首先，您应该给我们一个 AC 函数图。您显示的 PSD 的逆 FT 将是周期性音调的线性组合。

其次，尝试移除窗口，将其做成一个盒子，或者如果可以的话跳过这一步。

第三，尝试用 FFT 替换 DFT（我只浏览了 fftw3 库文档，也许这是一个选项）。

最后，尝试输入白噪声。您可以使用伯努利分布，也可以只使用高斯分布。 AC 将是一个 delta 函数，尽管样本 AC 不是。这应该会给你一个（样本）白色 PSD 分布。

我希望这些建议对您有所帮助。