c代码实现 ifft运算_X^n+1=0上的FFT和IFFT（基2）——C语言实现

我们一般意义上学习的FFT都是基于

的，即FFT中的单位根我们取的是

，但是在某些情况下我们需要

上的FFT和IFFT变换。

1、直接想到的思路是把

的根替换成

的根。

解法：

的根可以使用

的2n个根中的奇数次根得到，即

，但是这种做法在FFT运算中可行，在IFFT逆运算下则不可行，我们一般的IFFT运算时把

替换成

，并且最后除以一个n得到IFFT运算的结果。如下

但是我们需要在

上做IFFT变换的时候不能简单的把根替换成

，因为根据FFT的点值多项式的形式，只有根是

的形式的时候，才可以使用

因为根是

 的形式的时候，

在

上IFFT求逆的时候，

不成立，直接替换根的做法是不可行的。

2、新的做法，扩展到2n次

设

,

,

,f(x)是n次多项式。

令:

则有：

FFT计算步骤：

(1)计算

(2)此时F(x)是m次方的，计算F(x)在

上的FFT(就是以前一般形式的FFT)

(3)输出F(x) 的FFT变换之后的奇数项，即为f(x)在

上的FFT结果

IFFT计算步骤：

设

是f(x)在

上的FFT结果

(1)将

扩展，使其奇数项为

，偶数项为0，扩展到2n次

(2)使用2n阶IFFT求出扩展后多项式的逆变换的值

(3)设(2)中逆变换对应的扩展多项式逆变换为

,令

还原出来的n次

大致整体思路就是扩展到2n次，然后使用

上的FFT和IFFT求出

上的FFT和IFFT变换

下面贴C语言代码：

#include "pch.h"

#define _CRT_SECURE_NO_WARNINGS

#include "stdlib.h"

#include "math.h"

#include "stdio.h"

#define N 8

#define MAXN 100

#define Pi 3.1415927 //定义圆周率Pi

#define LEN sizeof(struct Compx) //定义复数结构体大小

//-----定义复数结构体-----------------------

typedef struct Compx

{

double real;

double imag;

}Compx;

//-----复数乘法运算函数---------------------

struct Compx mult(struct Compx b1, struct Compx b2)

{

struct Compx b3;

b3.real = b1.real*b2.real - b1.imag*b2.imag;

b3.imag = b1.real*b2.imag + b1.imag*b2.real;

return(b3);

}

//-----复数减法运算函数---------------------

struct Compx sub(struct Compx a, struct Compx b)

{

struct Compx c;

c.real = a.real - b.real;

c.imag = a.imag - b.imag;

return(c);

}

//-----复数加法运算函数---------------------

struct Compx add(struct Compx a, struct Compx b)

{

struct Compx c;

c.real = a.real + b.real;

c.imag = a.imag + b.imag;

return(c);

}

void fft(Compx *a, int n, int inv);

void fft_1(Compx *a, int n, int inv);

int main()

{

int i;

int x[N] = { 0 }, y[N] = { 0 };

printf("\nN=%d\n", N);

printf("\n输入两个多项式的系数，输入系数为N多项式长度的一半\n");

printf("\n输入第一个多项式的系数\n");

for (i = 0; i < N / 2; i++)

{

scanf("%d", &x[i]);

}

struct Compx * a = (struct Compx *)malloc(N*LEN);//为结构体分配存储空间

struct Compx * b = (struct Compx *)malloc(N*LEN);

struct Compx * c = (struct Compx *)malloc(N*LEN);

struct Compx * F = (struct Compx *)malloc(2*N*LEN);

//初始化=======================================

printf("\na多项式初始化：\n");

for (i = 0; i < N; i++)

{

a[i].real = x[i];

a[i].imag = 0;

printf("%.4f ", a[i].real);

printf("+%.4fj ", a[i].imag);

printf("\n");

}

/*--------------x^2n-1=0的解法----start----------*/

printf("\n--------------------------FFT---------------------------------\n");

int m = 2 * N;

int n = N;

double f[2 * N] = { 0 };

for (i = 0; i < N; i++) {

f[i] = 0.5 * x[i];

}

for (i = N; i < 2*N; i++) {

f[i] = -0.5 * x[i-N];

}

printf("\nf多项式初始化：\n");

for (i = 0; i < 2*N; i++)

{

F[i].real = f[i];

F[i].imag = 0;

printf("%.4f ", F[i].real);

printf("+%.4fj ", F[i].imag);

printf("\n");

}

fft(F, m, 1);

printf("\nF多项式FFT计算结果：\n");

for (i = 0; i < 2*N; i++)

{

printf("%.4f ", F[i].real);

printf("+%.4fj ", F[i].imag);

printf("\n");

}

for (i = 0; i < N; i++ ) {

a[i] = F[2 * i + 1];

}

printf("\n--------------------------IFFT---------------------------------\n");

fft(F, m, -1);

for (i = 0; i < m; i++) {

F[i].real = F[i].real / m;

F[i].imag = F[i].imag / m;

}

printf("\nF多项式IFFT计算结果：\n");

for (i = 0; i < 2 * N; i++)

{

printf("%.4f ", F[i].real);

printf("+%.4fj ", F[i].imag);

printf("\n");

}

//F(x)的低次

double temp_low[N] = { 0 };

double temp_high[N] = { 0 };

for (i = 0; i < N; i++)

{

temp_low[i] = F[i].real;

}

for (i = N; i < 2*N; i++)

{

temp_high[i-N] = F[i].real;

}

double res[N] = { 0 };

for (i = 0; i < N; i++)

{

res[i] = temp_low[i] - temp_high[i];

}

printf("\nIFFT计算结果：\n");

for (i = 0; i < N; i++)

{

printf("%.4f", res[i]);

printf("\n");

}

/*--------------x^2n+1=0的解法----end----------*/

//fft_1(a, N, 1);

printf("\n第一个多项式FFT计算结果：\n");

for (i = 0; i < N; i++)

{

printf("%.4f ", a[i].real);

printf("+%.4fj ", a[i].imag);

printf("\n");

}

return 0;

}

//x^n+1=0的FFT形式

void fft_1(Compx *a, int n, int inv) {

if (n == 1)return;

int mid = n / 2;

static Compx b[MAXN];

int i;

for (i = 0; i < mid; i++) {

b[i] = a[i * 2];

b[i + mid] = a[i * 2 + 1];

}

for (i = 0; i < n; i++)

a[i] = b[i];

fft(a, mid, inv);

fft(a + mid, mid, inv);//分治

for (i = 0; i < mid; i++)

{

Compx x;

x.real = cos(2 * Pi*(2*i+1) / (2*n));

x.imag = inv * sin(2 * Pi*(2 * i + 1) / (2*n));

b[i] = add(a[i], mult(x, a[i + mid]));

b[i + mid] = sub(a[i], mult(x, a[i + mid]));

}

for (i = 0; i < n; i++)

a[i] = b[i];

}

//x^n-1=0的FFT形式

void fft(Compx *a, int n, int inv) {

if (n == 1)return;

int mid = n / 2;

static Compx b[MAXN];

int i;

for (i = 0; i < mid; i++) {

b[i] = a[i * 2];

b[i + mid] = a[i * 2 + 1];

}

for (i = 0; i < n; i++)

a[i] = b[i];

fft(a, mid, inv);

fft(a + mid, mid, inv);//分治

for (i = 0; i < mid; i++)

{

Compx x;

x.real = cos(2 * Pi*i / n);

x.imag = inv * sin(2 * Pi*i / n);

b[i] = add(a[i], mult(x, a[i + mid]));

b[i + mid] = sub(a[i], mult(x, a[i + mid]));

}

for (i = 0; i < n; i++)

a[i] = b[i];

}