在matlab中用C语言设计FFT

在 MATLAB 中使用 C 语言编写 FFT 有两种方法：

1. 使用 MATLAB 的内置函数 coder，可以将 MATLAB 代码转换为 C 代码。

示例代码：

function y = my_fft(x)
%#codegen
coder.cinclude('fftw3.h');
coder.cinclude('fftw3_mkl.h');

N = length(x);
y = zeros(size(x));
p = coder.opaque('void *', 'NULL');
q = coder.opaque('void *', 'NULL');

if N > 1
    p = coder.opaque('fftw_plan', 'NULL');
    q = coder.opaque('fftw_plan', 'NULL');
    x = fftw_mkl(x);
    y = fftw_mkl(y);

    p = coder.ceval('fftw_plan_dft_1d', N, coder.ref(x), coder.ref(x), FFTW_FORWARD, FFTW_ESTIMATE);
    q = coder.ceval('fftw_plan_dft_1d', N, coder.ref(y), coder.ref(y), FFTW_FORWARD, FFTW_ESTIMATE);

    coder.ceval('fftw_execute_dft', p, coder.ref(x), coder.ref(x));
    coder.ceval('fftw_execute_dft', q, coder.ref(y), coder.ref(y));

    coder.ceval('fftw_destroy_plan', p);
    coder.ceval('fftw_destroy_plan', q);
    coder.ceval('fftw_cleanup');
end

2. 使用 MATLAB 的 MEX 文件。

示例代码：

#include "mex.h"
#include "fftw3.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *xreal, *ximag;
    fftw_complex *x;
    fftw_plan plan;
    int N;

    /* Check for proper number of arguments */
    if (nrhs != 1) {
        mexErrMsgIdAndTxt("MATLAB:fftw_mex:invalidNumInputs",
            "One input argument required.");
    } else if (nlhs > 1) {
        mexErrMsgIdAndTxt("MATLAB:fftw_mex:maxlhs",
            "Too many output arguments.");
    }

    /* Check type of input argument */
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
        mexErrMsgIdAndTxt("MATLAB:fftw_mex:inputNotReal",
            "Input argument must be real double array.");
    }

    /* Get size of input argument */
    N = mxGetNumberOfElements(prhs[0]);

    /* Allocate memory for input and output arrays */
    xreal = mxGetPr(prhs[0]);
    ximag = mxMalloc(N * sizeof(double));
    x = fftw_malloc(N * sizeof(fftw_complex));

    /* Fill complex array with input data */
    for (int i = 0; i < N; i++) {
        x[i][0] = xreal[i];
        x[i][1] = 0.0;
    }

    /* Create FFT plan */
    plan = fftw_plan_dft_1d(N, x, x, FFTW_FORWARD, FFTW_ESTIMATE);

    /* Execute FFT plan */
    fftw_execute(plan);

    /* Copy results to output array */
    plhs[0] = mxCreateDoubleMatrix(N, 1, mxCOMPLEX);
    xreal = mxGetPr(plhs[0]);
    ximag = mxGetPi(plhs[0]);
    for (int i = 0; i < N; i++) {
        xreal[i] = x[i][0];
        ximag[i] = x[i][1];
    }

    /* Clean up memory and FFT plan */
    fftw_destroy_plan(plan);
    fftw_free(x);
    mxFree(ximag);
}

这两种方法都需要使用 FFTW 库，所以要先安装 FFTW 库。