C++用代码验证“一切函数皆可傅里叶”

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


#define pi 3.1415926
#define rows 3
#define colums 5

typedef struct {
    float re;// really 
    float im;// imaginary 
}complex, * pcomplex;

complex complexadd(complex a, complex b) //复数加
{
    complex rt;
    rt.re = a.re + b.re;
    rt.im = a.im + b.im;
    return rt;
}

complex complexMult(complex a, complex b) //复数乘
{
    complex rt;
    rt.re = a.re * b.re - a.im * b.im;
    rt.im = a.im * b.re + a.re * b.im;
    return rt;
}
//离散傅里叶变换
void dft(complex F[], complex x[], int N) //X[]标识变换后频域，x[]为时域采样信号
{
    complex temp;
    int k, n;
    for (int k = 0; k < N; k++) {
        F[k].re = 0;
        F[k].im = 0;
        for (int n = 0; n < N; n++) {
            temp.re = (float)cos(2 * pi * k * n / N);
            temp.im = -(float)sin(2 * pi * k * n / N);
            F[k] = complexadd(F[k], complexMult(x[n], temp));
        }
        F[k].re = F[k].re / N;
        F[k].im = F[k].im / N;
    }
}

int main() {
    complex samples[colums], X[colums]; //samples[]示例
    for (int i = 0; i < colums; i++) {
        samples[i].re = i+1;
        samples[i].im = 0;
    }
    printf("原函数:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f,%f)\n", samples[i].re, samples[i].im);
    }
    printf("\n");

    dft(X, samples, colums);
    printf("原函数傅里叶变换 DFI:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f,%f)\n", X[i].re, X[i].im);
    }
    printf("\n");

    complex verify[colums];
    for (int i = 0; i < colums; i++) {
        verify[i].re = 0;
        verify[i].im = 0;
    }
    float Amplitude = 0;//幅度
    float Phase = 0;//相位
    float Frequency = 0;//频率
    float t = 0;//自变量
    float Re = 0;
    float Im = 0;
    printf("幅度      频率     采样弧度    相位    实部      虚部 \n");
    for (int j = 0; j < colums; j++){
        Frequency = j;//计算每个三角函数的频率
        for (int i = 0; i < colums; i++) {
            Amplitude = sqrt(pow(X[i].re, 2) + pow(X[i].im, 2));//计算每个三角函数的幅度
            Phase = atan2(X[i].im, X[i].re);//计算每个三角函数的相位
            t = 2 * pi * i / colums;//计算采样弧度
            Re = Amplitude * cos(Frequency * t + Phase);//计算函数的实部
            Im = Amplitude * sin(Frequency * t + Phase);//计算函数的虚部
            verify[j].re += Re;
            verify[j].im += Im;

            //verify[j].re += Amplitude * cos(2*pi * j * i /colums  + atan2(X[i].im, X[i].re));//整体计算误差会更小
            //verify[j].im += Amplitude * sin(2*pi * j * i /colums  + atan2(X[i].im, X[i].re));
            printf("(%f,%f,%f,%f,%f,%f)\n", Amplitude, Frequency,t, Phase, Re, Im);
        }
    }
    printf("\n");

    printf("验证结果:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f,%f)\n", verify[i].re, verify[i].im);
    }
    printf("\n");
}