python - fftw 中的真实一维 DFT

Question

我正在尝试在我的 C++ 代码中使用 FFTW3，并且我想使用 scipy.fftpack.fft 在 python 中完成同样的事情 用于真正的一维信号。我只是制作一维信号并从信号中找到频率。我不知道我哪里错了。 C++ 代码给我奇怪的结果。

这是scipy的例子:

    from scipy.fftpack import fft
    # Number of sample points
    N = 600
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N)
    y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
    yf = fft(y)
    xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
    import matplotlib.pyplot as plt
    plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
    plt.grid()
    plt.show()

and here is my c++ code:(**Edited**)

#include <fftw3.h>

#define NUM_POINTS 1024
/* Never mind this bit */
#include <cstdio>
#include <cmath>
#include <iostream>
#include <fstream>

#define REAL 0
#define IMAG 1

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

    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)*fabs(result[i][REAL]);
        f2<< (double)i << " "<<yf <<"\n";
    }
    f2.close();
    fftw_destroy_plan(plan);
    delete[] signal,theta;

    return 0;
}
======================================================
import pylab as pl 
import numpy as np 

signal = np.loadtxt("signal.txt")
pl.plot(signal[:,0],signal[:,1])

pl.figure()
result = np.loadtxt("result.txt")
pl.plot(result[:,0],result[:,1])

pl.show()

Answer 1

首先，您分配的内存太多:您只需要 fftw_complex result[NUM_POINTS/2+1]。但这不会破坏您的输出。

问题是fabs(result[i][REAL])，是result[i]实部的范数。在您的示例中，结果没有实部，只有虚部。函数 numpy.abs 采用复范数。