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

Question

我想用至少 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 中的频率？

Answer 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

...

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

calculate the precise, non integer frequencies from time series in (FFTW)

c++

fftw