C++ 通过频谱怎么计算功率谱

在++ 中，可以使用快速傅里叶变换（FFT）库（如 FFTW）来计算信号的频谱。具体计算方法如下：

假设有一个长度为 $N$ 的实数信号 $x[n]$（$n$ 为时间序列的索引）。首先对 $x[n]$ 进行 $N$ 点 FFT 变换得到其频谱 $X[k]$，其中 $k$ 为频率序列的索引。然后可以计算每个频率点 $f_k = k\Delta f$ 的功率谱密度 $S_k$，其中 $\Delta f = \frac{1}{N\Delta t}$，$\Delta t$ 为采样时间间隔。具体计算公式为：

$$S_k = \frac{1}{N} \cdot \frac{|X[k]|^2}{\Delta f}$$

在 C++ 中，可以使用 FFTW 库来计算 FFT 和功率谱。以下是一个示例代码，假设输入信号为一个名为 `input` 的数组，长度为 `N`，采样时间间隔为 `dt`：

#include <cmath>
#include <fftw3.h>

// 计算幅值平方
inline double square(double x) { return x * x; }

// 计算功率谱密度
inline double power_spectrum(double re, double im, double dt, int N) {
    double df = 1.0 / (N * dt);
    return square(re) + square(im) / N / df;
}

int main() {
    int N = 1024; // 信号长度
    double dt = 0.01; // 采样时间间隔

    // 分配输入和输出缓冲区
    double* input = new double[N];
    fftw_complex* output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));

    // 创建 FFT 句柄
    fftw_plan plan = fftw_plan_dft_r2c_1d(N, input, output, FFTW_ESTIMATE);

    // 执行 FFT 变换
    fftw_execute(plan);

    // 计算功率谱密度
    double* psd = new double[N / 2 + 1];
    for (int k = 0; k <= N / 2; k++) {
        psd[k] = power_spectrum(output[k][0], output[k][1], dt, N);
    }

    // 释放资源
    fftw_destroy_plan(plan);
    fftw_free(output);
    delete[] input;
    delete[] psd;

    return 0;
}

这段代码中，`fftw_plan_dft_r2c_1d` 函数创建了一个实数到复数的 FFT 句柄，`fftw_execute` 函数执行 FFT 变换，将结果存储在 `output` 缓冲区中。然后遍历 `output` 缓冲区，计算每个频率点的功率谱密度。注意到这里只需要计算前一半频率点，因为 FFT 的结果是对称的。最后释放创建的资源，包括 FFT 句柄、输出缓冲区和输入缓冲区。