#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

typedef struct {
    int    N;
    float *W;
    float *T;
    int   *order;
    int    ifft;
} FFT_CONTEXT;

static void complex_add(float *r, float *c1, float *c2) { r[0] = c1[0] + c2[0]; r[1] = c1[1] + c2[1]; }
static void complex_sub(float *r, float *c1, float *c2) { r[0] = c1[0] - c2[0]; r[1] = c1[1] - c2[1]; }
static void complex_mul(float *r, float *c1, float *c2) { r[0] = c1[0] * c2[0] - c1[1] * c2[1]; r[1] = c1[1] * c2[0] + c1[0] * c2[1]; }

static int reverse_bits(int n)
{
    n = ((n & 0xAAAAAAAA) >> 1 ) | ((n & 0x55555555) << 1 );
    n = ((n & 0xCCCCCCCC) >> 2 ) | ((n & 0x33333333) << 2 );
    n = ((n & 0xF0F0F0F0) >> 4 ) | ((n & 0x0F0F0F0F) << 4 );
    n = ((n & 0xFF00FF00) >> 8 ) | ((n & 0x00FF00FF) << 8 );
    n = ((n & 0xFFFF0000) >> 16) | ((n & 0x0000FFFF) << 16);
    return n;
}

static void fft_execute_internal(FFT_CONTEXT *ctxt, float *data, int n, int w)
{
    int  i;
    for (i=0; i<n/2; i++) {
        // C = (A + B)
        // D = (A - B) * W
        float  A[2] = { data[(0   + i) * 2 + 0], data[(0   + i) * 2 + 1] };
        float  B[2] = { data[(n/2 + i) * 2 + 0], data[(n/2 + i) * 2 + 1] };
        float *C    = &(data[(0   + i) * 2]);
        float *D    = &(data[(n/2 + i) * 2]);
        float *W    = &(ctxt->W[i*w*2]);
        float  T[2];
        complex_add(C, A, B);
        complex_sub(T, A, B);
        complex_mul(D, T, W);
    }

    n /= 2;
    w *= 2;
    if (n > 1) {
        fft_execute_internal(ctxt, data + 0    , n, w);
        fft_execute_internal(ctxt, data + n * 2, n, w);
    }
}

void fft_free(void *c)
{
    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c;
    if (!ctxt) return;
    if (ctxt->W    ) free(ctxt->W    );
    if (ctxt->T    ) free(ctxt->T    );
    if (ctxt->order) free(ctxt->order);
    free(ctxt);
}

void *fft_init(int n, int ifft)
{
    int shift, i;

    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)calloc(1, sizeof(FFT_CONTEXT));
    if (!ctxt) return NULL;

    ctxt->N     = n;
    ctxt->ifft  = ifft;
    ctxt->W     = (float*)calloc(n, sizeof(float) * 1);
    ctxt->T     = (float*)calloc(n, sizeof(float) * 2);
    ctxt->order = (int  *)calloc(n, sizeof(int  ) * 1);
    if (!ctxt->W || !ctxt->T || !ctxt->order) {
        fft_free(ctxt);
        return NULL;
    }

    for (i=0; i<ctxt->N/2; i++) {
        ctxt->W[i * 2 + 0] =(float) cos((ctxt->ifft ? -2 : 2) * M_PI * i / ctxt->N);
        ctxt->W[i * 2 + 1] =(float)-sin((ctxt->ifft ? -2 : 2) * M_PI * i / ctxt->N);
    }

    shift = 32 - (int)ceil(log(n)/log(2));
    for (i=0; i<ctxt->N; i++) {
        ctxt->order[i] = (unsigned)reverse_bits(i) >> shift;
    }
    return ctxt;
}

void fft_execute(void *c, float *in, float *out)
{
    int i;
    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c;
    memcpy(ctxt->T, in, sizeof(float) * 2 * ctxt->N);
    fft_execute_internal(ctxt, ctxt->T, ctxt->N, 1);
    for (i=0; i<ctxt->N; i++) {
        out[ctxt->order[i] * 2 + 0] = ctxt->T[i * 2 + 0] / (ctxt->ifft ? ctxt->N : 1);
        out[ctxt->order[i] * 2 + 1] = ctxt->T[i * 2 + 1] / (ctxt->ifft ? ctxt->N : 1);
    }
}

static void digitalstr2complex(float *complex, int len, char *str)
{
    int n = strlen(str), i;
    if (n < len) memset(complex + n * 2, 0, (len - n) * sizeof(float) * 2);
    for (i=n-1; i>=0; i--) {
        if (str[i] >= '0' && str[i] <= '9') complex[0] = str[i] - '0';
        else complex[0] = 0;
        complex[1] = 0;
        complex   += 2;
    }
}

int main(int argc, char *argv[])
{
    #define FFTLEN 512
    char *stra = "12341234234";
    char *strb = "234594350";
    void *fftf, *ffti;
    float fa[FFTLEN * 2], fb[FFTLEN * 2], fc[FFTLEN * 2];
    int   out[FFTLEN], carry, flag = 0, i;

    if (argc > 2) { stra = argv[1]; strb = argv[2]; }
    digitalstr2complex(fa, FFTLEN, stra);
    digitalstr2complex(fb, FFTLEN, strb);

    fftf = fft_init(FFTLEN, 0);
    ffti = fft_init(FFTLEN, 1);
    if (fftf && ffti) {
        fft_execute(fftf, fa, fa);
        fft_execute(fftf, fb, fb);
        for (i=0; i<FFTLEN; i++) complex_mul(fc + i * 2, fa + i * 2, fb + i * 2);
        fft_execute(ffti, fc, fc);
    }
    fft_free(fftf);
    fft_free(ffti);

    for (carry=0,i=0; i<FFTLEN; i++) {
        out[i] = (int)(fc[i * 2] + 0.5) + carry;
        carry  = out[i] / 10;
        out[i]%= 10;
    }
    for (i=FFTLEN-1; i>=0; i--) {
        if (flag == 0) flag = !!out[i];
        if (flag) printf("%d", out[i]);
    }
    printf("\n");
    return 0;
}

[此贴子已经被作者于2021-7-23 18:19编辑过]

这个算法的原理很简单，卷积定理：时域的卷积等价于频域的乘法
这样可以把 O(n * n) 的复杂度简化为 O(n * log(n)) 的复杂度
代码也就 100 多行就搞定了，简单高效