从 (FFTW) 中的时间序列计算精确的非整数频率答案

【问题标题】：calculate the precise, non integer frequencies from time series in (FFTW)从 (FFTW) 中的时间序列计算精确的非整数频率
【发布时间】：2018-02-19 10:33:02
【问题描述】：

我想用至少 3 个十进制值精确计算时间序列的频率。这是一个计算整数值频率的简单示例。

#include <fftw3.h>
#include <cstdio>
#include <cmath>
#include <iostream>
#include <fstream>

#define REAL 0
#define IMAG 1
#define NUM_POINTS 1024


void acquire_signal(double *signal, double *theta) {
    /* Generate two sine waves of different frequencies and
     * amplitudes.
     */

    int i;
    for (i = 0; i < NUM_POINTS; ++i) {
        theta[i] = (double)i / (double)NUM_POINTS;
        signal[i] = 1.0*sin(50.0 * 2.0 * M_PI * theta[i]) +
                    0.5*sin(80.0 * 2.0 * M_PI * theta[i]);
    }
}

int main() {
    unsigned flags{0};
    double *theta  = new double[NUM_POINTS];
    double *signal = new double[NUM_POINTS];

    fftw_complex result[NUM_POINTS/2+1];

    fftw_plan plan = fftw_plan_dft_r2c_1d(NUM_POINTS,
                                         signal,
                                         result,
                                         flags);
    acquire_signal(signal,theta);
    fftw_execute(plan);

    //save signal and result
    std::ofstream f1,f2;
    f1.open ("signal.txt");
    for (int i=0; i<NUM_POINTS; i++){
        f1 <<theta[i]<<" "<<signal[i]<<"\n";
    }

    f1.close();
    f2.open("result.txt");

    for (int i=0; i<NUM_POINTS/2; i++){
        double yf = 2.0/(double)(NUM_POINTS)* sqrt(result[i][REAL]*result[i][REAL]+ result[i][IMAG]* result[i][IMAG]);
        f2<< (double)i << " "<<yf <<"\n";
    }
    f2.close();
    fftw_destroy_plan(plan);
    delete[] signal,theta;

    return 0;
}

但是如果我有，我应该如何更改代码

signal = 1.0*sin(50.350 * 2.0 * M_PI * theta[i]) +
         0.5*sin(80.455 * 2.0 * M_PI * theta[i]);

改变时间和频率的单位是否合适？例如1000*s 中的时间和kHz 中的频率？

【问题讨论】：

这个问题可能是关于数学而不是 c++ 或编码的吗？
是的，我应该把问题移到右边？
您的频率取决于您的采样频率和样本数量。您不能只调整精度。您需要更多点或不同的采样频率
只需更改数字即可将您的线路从 50 和 80 赫兹转移到 50.350 和 80.455 赫兹，假设您有 1024 行乘以 1024 赫兹。但是你仍然有 1Hz 的分辨率。
无论时间序列越长，我最终都会用f2<< (double)i << " "<<yf; 打印频率和幅度，这是一个整数，表示频率。 i 在 [0,NUM_POINTS/2] 中。我错了吗？

标签： c++ fftw

【解决方案1】：

只需更改 sin 中的数字，您的线路将从 50 和 80 转移到 50.350 和 80.455 Hz，并假设您有 1024 行乘以 1024 Hz。但是你仍然有 1Hz 的分辨率。在相同的采样频率下，您需要更多行 (x1000) 才能获得更大的分辨率。

例如，如果您想要 1/4 Hz 分辨率，则需要 4 倍以上的样本，因此按照 1024 Hz 采样率，您需要 fs * 4 样本：

...

#define NUM_POINTS (1024 * 4)
double fs = 1024;  // Sample rate in Hz

void acquire_signal(double *signal, double *theta) {
    /* Generate two sine waves of different frequencies and
     * amplitudes.
     */

    int i;
    for (i = 0; i < NUM_POINTS; ++i) {
        theta[i] = (double)i / (double)fs;
        signal[i] = 1.0*sin(50.0 * 2.0 * M_PI * theta[i]) +
                    0.5*sin(80.0 * 2.0 * M_PI * theta[i]);
    }
}
....
for (int i=0; i< (NUM_POINTS/2 + 1) ; i++){
    double yf = 2.0/(double)(NUM_POINTS)* sqrt(result[i][REAL]*result[i][REAL]+ result[i][IMAG]* result[i][IMAG]);
    f2 << (double)i * fs / ( NUM_POINTS ) << " "<<yf <<"\n";
}

0 2.90715e-16

0.25 1.19539e-16

0.5 2.15565e-16

0.75 2.88629e-16

1 3.05084e-16

1.25 3.864e-16

...

49.75 9.47968e-16

50 1

50.25 1.12861e-15

50.5 4.95946e-16

50.75 6.9016e-16

...

【讨论】：