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));
        }
    }
}
//离散傅里叶逆变换
void idft(complex X[], complex x[], int N)
{
    complex temp;
    int k, n;
    for (int k = 0; k < N; k++)
    {
        x[k].re = 0;
        x[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);
            x[k] = complexadd(x[k], complexMult(X[n], temp));
        }
        x[k].re /= N;
        x[k].im /= N;
    }
}

void dft2d(complex F[][colums], complex x[][colums], int R, int C) //二维离散傅里叶变换
{
    complex temp;
    int k, n, i, j;
    for (i = 0; i < R; i++)
    {
        for (j = 0; j < C; j++)
        {
            F[i][j].re = 0;
            F[i][j].im = 0;
            for (k = 0; k < R; k++)
            {
                for (n = 0; n < C; n++)
                {
                    temp.re = (float)cos(2 * pi * (i * k / (float)R + j * n / (float)C));
                    temp.im = -(float)sin(2 * pi * (i * k / (float)R + j * n / (float)C));
                    F[i][j] = complexadd(F[i][j], complexMult(x[k][n], temp));
                }
            }
        }
    }
}

void idft2d(complex X[][colums], complex x[][colums], int R, int C) //二维离散傅里叶逆变换
{
    complex temp;
    int k, n, i, j;
    for (i = 0; i < R; i++)
    {
        for (j = 0; j < C; j++)
        {
            x[i][j].re = 0;
            x[i][j].im = 0;
            for (k = 0; k < R; k++)
            {
                for (n = 0; n < C; n++)
                {
                    temp.re = (float)cos(2 * pi * (i * k / (float)R + j * n / (float)C));
                    temp.im = (float)sin(2 * pi * (i * k / (float)R + j * n / (float)C));
                    x[i][j] = complexadd(x[i][j], complexMult(X[k][n], temp));
                }
            }
            x[i][j].re /= R * C;
            x[i][j].im /= R * C;
        }
    }
}
int main()
{
    complex samples[colums], X[colums], x[colums]; //samples[]示例

    for (int i = 0; i < colums; i++)
    {
        samples[i].re = i;
        samples[i].im = 0;
    }
    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");

    idft(X, x, colums);
    printf("一维傅里叶逆变换 IDFI:\n");
    for (int i = 0; i < colums; i++)
    {
        printf("(%f,%f)\n", x[i].re, x[i].im);
    }
    printf("\n");

    //二维数组的傅里叶变换示例
    complex samples2d[rows][colums], X2d[rows][colums], x2d[rows][colums]; //samples2d[]示例

    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < colums; j++)
        {
            samples2d[i][j].re = i + j;
            samples2d[i][j].im = 0;
        }
    }

    dft2d(X2d, samples2d, rows, colums);
    printf("二维傅里叶变换 2D DFI:\n");
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < colums; j++)
        {
            printf("(%f,%f)\t", X2d[i][j].re, X2d[i][j].im);
        }
        printf("\n");
    }
    printf("\n");

    idft2d(X2d, x2d, rows, colums);
    printf("二维傅里叶逆变换  2D IDFI:\n");
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < colums; j++)
        {
            printf("(%f,%f)\t", x2d[i][j].re, x2d[i][j].im);
        }
        printf("\n");
    }
}