# 快速傅里叶变换

## 准备

FT (Fourier Transform): 傅里叶变换

DFT (Discrete Fourier Transform): 离散傅里叶变换

FFT (Fast Fourier Transform): 快速傅里叶变换

IFFT (Inverse Fourier Transform): 快速傅里叶变换逆变换

### 多项式的点值表示

$$\left[\begin{matrix}A(x_0)\\A(x_1)\\ \vdots\\A(x_n)\end{matrix}\right]=\left[\begin{matrix}1&x_0&x_0^2&\cdots&x_0^{n}\\1&x_1&x_1^2&\cdots&x_1^n\\&&\vdots&&\\1&x_n&x_n^2&\cdots&x_n^n\end{matrix}\right]\left[\begin{matrix}a_0\\a_1\\\vdots \\a_n\end{matrix}\right]$$

### 单位根

$$\omega_n^k=(\cos{\frac{k}{n}2\pi},\sin{\frac{k}{n}2\pi})$$

\begin{aligned}\omega_n^k\times\omega_n^1=&(\cos{\frac{k}{n}2\pi},\sin{\frac{k}{n}2\pi})*(\cos{\frac{1}{n}2\pi},\sin{\frac{1}{n}2\pi})\\=&(\cos{\frac{k}{n}2\pi}*\cos{\frac{1}{n}2\pi}-\sin{\frac{k}{n}2\pi}*\sin{\frac{1}{n}2\pi},\sin{\frac{k}{n}2\pi}*\cos{\frac{1}{n}2\pi}+\cos{\frac{k}{n}2\pi}*\sin{\frac{1}{n}2\pi})\\=&(\cos{\frac{k+1}{n}2\pi},\sin{\frac{k+1}{n}2\pi})\\=&\omega_n^{k+1}\end{aligned}

\begin{aligned}\omega_{n*x}^{k*x}=&(\cos{\frac{k*x}{n*x}2\pi},\sin{\frac{k*x}{n*x}2\pi})\\=&(\cos{\frac{k}{n}2\pi},\sin{\frac{k}{n}2\pi})\\=&\omega_n^k\end{aligned}

\begin{aligned}\omega_n^{k+\frac{n}{2}}=&(\cos{\frac{k+\frac{n}{2}}{n}}2\pi,\sin{\frac{k+\frac{n}{2}}{n}2\pi})\\=&(-\cos{\frac{k}{n}2\pi},-\sin{\frac{k}{n}2\pi})\\=&-(cos{\frac{k}{n}2\pi},sin{\frac{k}{n}2\pi})\\=&-\omega_n^k\end{aligned}
$$\omega_{\frac{n}{m}}^{k}=\omega_n^{mk}$$

\sum_{k=0}^{n-1}(\omega_n^{j-i})^k=\left\{\begin{aligned}0,\;i\neq j\\n,\;i=j\end{aligned}\right.

$A(x)=a_0x^0+a_1x^1+\cdots+a_{n-1}x^{n-1}$将$\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$作为$x$分别带入求值得到$(y_0,y_1,\cdots,y_{n-1})$, 将这些$y$作为系数,产生一个新的多项式$B(x)=b_0x^0+b_1x^1+\cdots+b_{n-1}x^{n-1}$, 将$\omega_n^{-0},\omega_n^{-1},\cdots,\omega_n^{-(n-1)}$作为$x$分别带入求值得到的$(z_0,z_1,\cdots,z_{n-1})$对于每个$z_k$, 有$z_k=a_k*n$

\begin{aligned}z_k=&\sum_{i=0}^{n-1}y_i*(\omega_n^{-k})^i\\=&\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j*(\omega_n^i)^j)*(\omega_n^{-k})^i\\=&\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j*(\omega_n^i)^j*(\omega^{-k})^i\\=&\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j*\omega_n^{i*j}*\omega_n^{-k*i}\\=&\sum_{j=0}^{n-1}a_j\left(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\right)\end{aligned}

1. 如果$j-k=0$那么$\omega_n^{j-k}$就是$\omega_n^0=1$, 所以$n$个$1$结果为$n$
2. 如果$j-k\neq 0$那么通过等比数列求和$(x^0+x^1+\cdots,x^{n-1}=\frac{x^n-1}{x-1})$可以发现, 结果$\frac{(\omega_n^{j-k})^n-1}{\omega_n^{j-k}-1}=\frac{(\omega_n^n)^{j-k}-1}{\omega_n^{j-k}-1}=0$

## FFT

### 蝶形运算

1. 左端线上有数字表示(没有则表示C和D为1):

1. 右端线上有数字表示(没有则表示C和D为1):

\begin{aligned}A(x)&=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+\cdots+a_{n-1}x^{n-1})\\&=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})\end{aligned}

\begin{aligned}A_0(x)&=a_0+a_2x^1+\cdots+a_{n-2}x^{\frac{n}{2}-1}\\A_1(x)&=a_1+a_3x^1+\cdots+a_{n-1}x^{\frac{n}{2}-1}\end{aligned}

\begin{aligned}A(\omega_n^m)&=A_0((\omega_n^m)^2)+\omega_n^mA_1((\omega_n^m)^2)\\&=A_0(\omega_n^{2m})+\omega_n^mA_1(\omega_n^{2m})\\&=A_{0}(\omega_{\frac{n}{2}}^m)+\omega_n^mA_1(\omega_{\frac{n}{2}}^m)\end{aligned}

\begin{aligned}A(\omega_n^{m+\frac{n}{2}})&=A_0(\omega_n^{2m+n})+\omega_n^{m+\frac{n}{2}} A_1(\omega_n^{2m+n})\\&=A_0(\omega_n^{2m}\omega_n^n)-\omega_n^mA_1(\omega_n^{2m}\omega_n^n)\\&=A_0(\omega_\frac{n}{2}^m)-\omega_n^mA_1(\omega_\frac{n}{2}^{m})\end{aligned}

\begin{aligned}A(\omega_n^m)&=A_0(\omega_\frac{n}{2}^m)+\omega_n^mA_1(\omega_\frac{n}{2}^m)\\A(\omega_n^{m+\frac{n}{2}})&=A_0(\omega_{\frac{n}{2}}^m)-\omega_n^mA_1(\omega_\frac{n}{2}^m)\end{aligned}

\begin{aligned}A(\omega^j)&=A_0(\omega^{2j})+\omega^jA_1(\omega^{2j})\\A(-\omega^j)&=A_0(\omega^{2j})-\omega^{j}A_1(\omega^{2j})\\j&\in\left\{0,1,\cdots,(n/2-1)\right\}\end{aligned}

### code

#define _USE_MATH_DEFINES
#include <bits/stdc++.h>

using namespace std;

int Extend(int num_coef);
void ReverseOrder(int n, float* coefR, float* coefI);
void Input(int num_coef, int n, float* dataR, float* dataI);
void Calculate(int P, int n, int r, int B, float* dataR, float* dataI);
void Print(int n, float* dataR, float* dataI);

int main(int argc, char** argv) {
cout.setf(ios_base::fixed, ios_base::floatfield);
cout.setf(ios_base::showpoint);
cout.precision(5);

int num_coef;
cin >> num_coef;
float *dataR, *dataI;
int n = Extend(num_coef);

dataR = new float[n];
dataI = new float[n];

Input(num_coef, n, dataR, dataI);

ReverseOrder(n, dataR, dataI);
for(int i=0;i<n;i++){
cout<<dataR[i]<<" ";
}
cout<<endl;

int M = log2(n);
for (int L = 1; L <= M; L++) {
int B = pow(2, L - 1);
for (int j = 0; j <= B-1; j++) {
int k = pow(2, M - L);
int P = j * k;
for (int i = 0; i <= k-1; i++) {
int r = j + 2 * B * i;
Calculate(P, n, r, B, dataR, dataI);
}
}
}
Print(n, dataR, dataI);
system("pause");
return 0;
}

int Extend(int num_coef) {
int m = log2(num_coef);
if (num_coef == pow(2, m))
return num_coef;
else
return pow(2, m + 1);
}

void Input(int num_coef, int n, float* dataR, float* dataI) {
cout << "R" << endl;
for (int i = 0; i < num_coef; i++) {
cin >> dataR[i];
}
cout << "I" << endl;
for (int i = 0; i < num_coef; i++) {
cin >> dataI[i];
}
if (n > num_coef)
for (int i = num_coef; i < n; i++) {
dataR[i] = 0;
dataI[i] = 0;
}
}

void ReverseOrder(int n, float* coefR, float* coefI) {
int m = log2(n);
int a1, an;
float tempR, tempI;
for (int i = 0; i < pow(2, m); i++) {
int j=0;
for (int k = 0; k < (m + 1) / 2; k++) {
a1 = 1;
an = pow(2, m - 1);
a1 <<= k;
an >>= k;
rear = i & a1;
j = a1|j;
if (0 != rear)
j = an|j;
}
if (i < j) {
swap(coefR[i], coefR[j]);
swap(coefI[i], coefI[j]);
}
}
}

void Calculate(int P, int n, int r, int B, float* dataR, float* dataI){
/*进行蝶形运算*/
//t = W * A1, W = cos(2 * PI * P / n) - i * sin(2 * PI * P / n)
cout<<M_PI<<endl;
float tR = dataR[r + B] * cos(2 * M_PI * P / n) + dataI[r + B] * sin(2 * M_PI * P / n);
float tI = dataI[r + B] * cos(2 * M_PI * P / n) - dataR[r + B] * sin(2 * M_PI * P / n);
//A(r) = A0 + W * A1
//A(r + B) = A0 - W * A1
dataR[r + B] = dataR[r] - tR;
dataI[r + B] = dataI[r] - tI;
dataR[r] = dataR[r] + tR;
dataI[r] = dataI[r] + tI;
}

void Print(int n, float* dataR, float* dataI){
cout << "result:" << endl;
for (int i = 0; i < n; i++) {
if(dataI[i] >= 0)
cout << dataR[i] << "+" << dataI[i] << "i" << "\t";
else
cout << dataR[i] << dataI[i] << "i" << "\t";
}
cout << endl;
}

## 并行FFT

1. 在每个进程中用mdataR[]mdataI[]来存储结果的实部和虚部
2. 在每级蝶形计算前先把上一级的结果广播给每个进程(从其他进程获取数据可能发生死锁)
3. 为被分配的元素由0号进程处理

1. 若对应鲽形组的上半部分，且鲽形组的下半部分也在此进程中，则直接更新两个元素：

2. 若对应鲽形组的上半部分，且蝶形组的下半部分不在此进程中，对下半部分的元素进行分析：

• 如果是已分配元素：更新当前元素，同时将下半部分对应元素发送给指定进程
• 如果是未分配的元素：更新当前元素，但是不发送下半部分对应的元素(发送给0号进程可能造成死锁)

注意：这里的发送的数据指的是计算的结果，并非蝶形组左边的数据，还望不要被图片所误导，这样标注仅是为了说明进程间的发送情况。

3. 若是对应蝶形组的下半部分，则接收其他进程发送过来的数据：

1. 若对应蝶形组的上半部分，则直接更新两个元素：
2. 若对应蝶形组的下半部分，则更新此元素：

### code

#define _USE_MATH_DEFINES
#include <math.h>
#include <mpi.h>
#include <stdlib.h>
#include <algorithm>
#include <iostream>

using namespace std;

int Extend(int num_coef);
void ReverseOrder(int n, float* coefR, float* coefI);
void Input(int num_coef, int n, float* dataR, float* dataI);
void Calculate(int P, int n, int r, int B, float* dataR, float* dataI, float& tempR, float& tempI, float& mdataR, float& mdataI);
void Print(int n, float* dataR, float* dataI);

int main(int argc, char** argv) {
cout.setf(ios_base::fixed, ios_base::floatfield);
cout.setf(ios_base::showpoint);
cout.precision(2);

int num_coef = atoi(argv[1]);
float *dataR, *dataI, *mdataR, *mdataI;
float tempR, tempI;

int myrank, size;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Status status;

int n = Extend(num_coef);
int msize = n / size;
mdataR = new float[msize];
mdataI = new float[msize];
dataR = new float[n];
dataI = new float[n];

if (0 == myrank) {
Input(num_coef, n, dataR, dataI);
ReverseOrder(n, dataR, dataI);
}
MPI_Barrier(MPI_COMM_WORLD);

int M = log2(n);
int start = myrank * msize;
for (int L = 1; L <= M; L++) {
MPI_Bcast(dataR, n, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(dataI, n, MPI_FLOAT, 0, MPI_COMM_WORLD);

int B = pow(2, L - 1);
int k = pow(2, M - L);

for (int j = 0; j < msize; j++) {
int r = start + j;
if ((r % (2 * B)) >= B) {
if (j - B < 0) {
MPI_Recv(&mdataR[j], 1, MPI_FLOAT, (r - B) / msize, r, MPI_COMM_WORLD, &status);
MPI_Recv(&mdataI[j], 1, MPI_FLOAT, (r - B) / msize, r, MPI_COMM_WORLD, &status);
}
continue;
}

int P = (r % B) * k;

Calculate(P, n, r, B, dataR, dataI, tempR, tempI, mdataR[j], mdataI[j]);
if (j + B + 1 > msize) {
if (r + B < msize * size) {
MPI_Send(&tempR, 1, MPI_FLOAT, (r + B) / msize, r + B, MPI_COMM_WORLD);
MPI_Send(&tempI, 1, MPI_FLOAT, (r + B) / msize, r + B, MPI_COMM_WORLD);
}
} else {
mdataR[j + B] = tempR;
mdataI[j + B] = tempI;
}
}

int remainStartID = msize * size;
if ((0 == myrank) && (remainStartID != n)) {
for (int j = 0; j < n - remainStartID; j++) {
int r = remainStartID + j;
int P = (r % B) * k;
if((r%(2*B))>=B){
if(j-B<0){
float tR = dataR[r] * cos(2 * M_PI * P / n) + dataI[r] * sin(2 * M_PI * P / n);
float tI = dataI[r] * cos(2 * M_PI * P / n) - dataR[r] * sin(2 * M_PI * P / n);
dataR[r] = dataR[r - B] - tR;
dataI[r] = dataI[r - B] - tI;
}
continue;
}
Calculate(P,n,r,B,dataR,dataI,dataR[r+B],dataI[r+B],dataR[r],dataI[r]);
}
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Gather(mdataR,msize,MPI_FLOAT,dataR,msize,MPI_FLOAT,0,MPI_COMM_WORLD);
MPI_Gather(mdataI,msize,MPI_FLOAT,dataI,msize,MPI_FLOAT,0,MPI_COMM_WORLD);

}
if(0==myrank){
Print(n,dataR,dataI);
}
MPI_Finalize();
return 0;
}

int Extend(int num_coef) {
int m = log2(num_coef);
if (num_coef == pow(2, m))
return (num_coef);
else
return (pow(2, m + 1));
}

void ReverseOrder(int n, float* coefR, float* coefI) {
int m = log2(n);
int a1, an;
for (int i = 0; i < pow(2, m); i++) {
int j = 0;
for (int k = 0; k < (m + 1) / 2; k++) {
a1 = 1;
an = pow(2, m - 1);
a1 <<= k;
an >>= k;
rear = i & a1;
j = a1 | j;
if (0 != rear)
j = an | j;
}
if (i < j) {
swap(coefR[i], coefR[j]);
swap(coefI[i], coefI[j]);
}
}
}

void Input(int num_coef, int n, float* dataR, float* dataI){
cout << "R: ";
int i;
for (int i = 0; i < num_coef; i++) {
cin >> dataR[i];
}
cout << "I: ";
for (i = 0; i < num_coef; i++) {
cin >> dataI[i];
}
if (n > num_coef)
for(; i < n; i++){
dataR[i] = 0;
dataI[i] = 0;
}
}

void Calculate(int P, int n, int r, int B, float* dataR, float* dataI, float& tempR, float& tempI, float& mdataR, float& mdataI) {
/*进行蝶形运算*/
//t = W * A1, W = cos(2 * PI * P / n) - i * sin(2 * PI * P / n)
float tR = dataR[r + B] * cos(2 * M_PI * P / n) + dataI[r + B] * sin(2 * M_PI * P / n);
float tI = dataI[r + B] * cos(2 * M_PI * P / n) - dataR[r + B] * sin(2 * M_PI * P / n);
//A(r) = A0 + W * A1
//A(r + B) = A0 - W * A1
tempR = dataR[r] - tR;
tempI = dataI[r] - tI;
mdataR = dataR[r] + tR;
mdataI = dataI[r] + tI;
}

void Print(int n, float* dataR, float* dataI) {
cout << "result:" << endl;
for (int i = 0; i < n; i++) {
if (dataI[i] >= 0)
cout << dataR[i] << "+" << dataI[i] << "i"
<< "\t";
else
cout << dataR[i] << dataI[i] << "i"
<< "\t";
}
cout << endl;
}