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


#define pi 3.1415926

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[]为时域采样信号，下同
{
    float data_re[10];
    float data_im[10];
    complex temp;
    complex data;
    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));
            data = complexadd(F[k], complexMult(x[n], temp));
            data_re[n] = temp.re;
            data_im[n] = temp.im;
        }
    }
}
//离散傅里叶逆变换
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[][5], complex x[][5], int N)
{
    complex temp;
    complex data;
    int k, n, i, j;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            F[i][j].re = 0;
            F[i][j].im = 0;
            for (k = 0; k < N; k++)
            {
                for (n = 0; n < N; n++)
                {
                    temp.re = (float)cos(2 * pi * (k * i + n * j) / N);
                    temp.im = -(float)sin(2 * pi * (k * i + n * j) / N);
                    F[i][j] = complexadd(F[i][j], complexMult(x[k][n], temp));
                    data = complexadd(F[i][j], complexMult(x[k][n], temp));
                }
            }
        }
    }
}
//二维数组的傅里叶逆变换
void idft2d(complex X[][5], complex x[][5], int N)
{
    complex temp;
    int k, n, i, j;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            x[i][j].re = 0;
            x[i][j].im = 0;
            for (k = 0; k < N; k++)
            {
                for (n = 0; n < N; n++)
                {
                    temp.re = (float)cos(2 * pi * (k * i + n * j) / N);
                    temp.im = (float)sin(2 * pi * (k * i + n * j) / N);
                    x[i][j] = complexadd(x[i][j], complexMult(X[k][n], temp));
                }
            }
            x[i][j].re /= N;
            x[i][j].im /= N;
        }
    }
}


#define N 5

int main()
{
    complex samples[N], X[N], x[N]; //samples[]示例

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

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

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

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

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

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