傅里叶变换
系数表示法与点值表示法
f ( x ) = a 0 ∗ x 0 + a 1 ∗ x 1 + a 2 ∗ x 2 + a 3 ∗ x 3 + . . . + a n ∗ x n f(x) = a_0 * x ^ 0 + a_1 * x^1 + a_2 * x^2 + a_3 * x^3+...+a_n * x ^ n f(x)=a0∗x0+a1∗x1+a2∗x2+a3∗x3+...+an∗xn
系数表示法：
(
a
0
,
a
1
,
a
2
,
a
3
,
a
4
,
.
.
.
,
a
n
)
(a_0,a_1,a_2,a_3,a_4,...,a_n)
(a0,a1,a2,a3,a4,...,an)
例如 ( 1 , 2 , 1 , 4 ) = > f ( x ) = 1 + 2 ∗ x 1 + 1 ∗ x 2 + 4 ∗ x 3 例如 (1,2,1,4) => f(x) = 1 + 2 * x^1 + 1*x^2 + 4 * x^3 例如(1,2,1,4)=>f(x)=1+2∗x1+1∗x2+4∗x3
点值表示法：n + 1个点坐标确定一条曲线
(
x
0
,
y
0
)
、
(
x
1
,
y
1
)
、
(
x
2
,
y
2
)
、
(
x
3
,
y
3
)
.
.
.
.
.
、
(
x
n
,
y
n
)
(x_0,y_0)、(x_1,y_1)、(x_2,y_2)、(x_3,y_3).....、(x_n,y_n)
(x0,y0)、(x1,y1)、(x2,y2)、(x3,y3).....、(xn,yn)
FFT算法的作用
f ( x ) = a 0 + a 1 ∗ x n 1 + a 2 ∗ x n 2 + a 3 ∗ x n 3 + ⋅ ⋅ ⋅ + a m ∗ x n m ( n = 1 , 2 , 3.... m ) ∑ 0 n f ( x ) f(x) = a_0 + a_1*x^1_n + a_2*x^2_n + a_3*x^3_n + ···+a_m*x^m_n (n = {1,2,3....m})\\ \sum^{n}_{0}{f(x)} f(x)=a0+a1∗xn1+a2∗xn2+a3∗xn3+⋅⋅⋅+am∗xnm(n=1,2,3....m)0∑nf(x)
等价于：
[ y 0 y 1 y 2 y 3 . . . y n ] = [ 1 x 0 x 0 2 x 0 3 ⋯ ⋯ , x 0 m 1 x 1 x 1 2 x 1 3 ⋯ ⋯ , x 1 m 1 x 2 x 2 2 x 2 3 ⋯ ⋯ , x 2 m 1 x 3 x 3 2 x 3 3 ⋯ ⋯ , x 3 m . . . 1 x n x n 2 x n 3 ⋯ ⋯ , x n m ] × [ a 0 a 1 a 2 a 3 . . . a n ] \left[ \begin{matrix} y_0\\ y_1\\ y_2\\ y_3\\ .\\ .\\ .\\ y_n \end{matrix} \right]= \left[ \begin{matrix} 1 &x_0 &x_0^2 &x_0^3 &\cdots &\cdots,x_0^m\\ 1 &x_1 &x_1^2 &x_1^3 &\cdots &\cdots,x_1^m\\ 1 &x_2 &x_2^2 &x_2^3 &\cdots &\cdots,x_2^m\\ 1 &x_3 &x_3^2 &x_3^3 &\cdots &\cdots,x_3^m\\ &&&.\\ &&&.\\ &&&.\\ 1 &x_n &x_n^2 &x_n^3 &\cdots &\cdots,x_n^m \end{matrix} \right]\times \left[ \begin{matrix} a_0\\ a_1\\ a_2\\ a_3\\ .\\ .\\ .\\ a_n \end{matrix} \right] y0y1y2y3...yn = 11111x0x1x2x3xnx02x12x22x32xn2x03x13x23x33...xn3⋯⋯⋯⋯⋯⋯,x0m⋯,x1m⋯,x2m⋯,x3m⋯,xnm × a0a1a2a3...an
FFT算法就是快速计算这个式子，由n个系数得到n个值构成点值表示
天才想法1-快速取值
思考：怎么在一个二项式里面得到一个六个点的值
解释：计算三次得到三个点，对称取值得到六个点
f ( x ) = a 0 ∗ x 0 + a 1 ∗ x 1 + a 2 ∗ x 2 + a 3 ∗ x 3 + a 4 ∗ x 4 + a 5 ∗ x 5 f(x) = a_0 * x ^ 0 + a_1 * x^1 + a_2 * x^2 + a_3 * x^3+a_4 * x ^ 4 +a_5 * x ^ 5 f(x)=a0∗x0+a1∗x1+a2∗x2+a3∗x3+a4∗x4+a5∗x5
思考：在一个多项式里面是否存在对称性，如果存在就只要求n / 2个点
解释：
第一、奇偶划分
f ( x ) = ( a 0 ∗ x 0 + a 2 ∗ x 2 + a 4 ∗ x 4 ) + ( a 1 ∗ x 1 + a 3 ∗ x 3 + a 5 ∗ x 5 ) f(x) = (a_0 * x ^ 0 + a_2 * x^2 +a_4 * x ^ 4) + (a_1 * x^1+ a_3 * x^3 +a_5 * x ^ 5)\\ f(x)=(a0∗x0+a2∗x2+a4∗x4)+(a1∗x1+a3∗x3+a5∗x5)
第二、多项式构造
P
0
(
y
=
x
2
)
=
(
a
0
+
a
2
∗
y
+
a
4
∗
y
)
P
1
(
y
=
x
2
)
=
(
a
1
+
a
3
∗
y
+
a
5
∗
y
2
)
f
(
x
)
=
P
0
(
x
2
)
+
x
P
1
(
x
2
)
得：
f
(
−
x
)
=
P
0
(
x
2
)
−
x
P
1
(
x
2
)
P_0(y=x^2)=(a_0 + a_2 * y + a_4 * y) \\ P_1(y=x^2)=(a_1 + a_3*y + a_5*y^2) \\ f(x) = P_0(x^2) + xP_1(x ^ 2)\\ 得：f(-x) = P_0(x^2) - xP_1(x ^ 2)
P0(y=x2)=(a0+a2∗y+a4∗y)P1(y=x2)=(a1+a3∗y+a5∗y2)f(x)=P0(x2)+xP1(x2)得：f(−x)=P0(x2)−xP1(x2)
第三、抽象成算法
f
(
x
)
=
P
0
(
x
2
)
+
x
P
1
(
x
2
)
f
(
−
x
)
=
P
0
(
x
2
)
−
x
P
1
(
x
2
)
f(x) = P_0(x^2) + xP_1(x ^ 2)\\ f(-x) = P_0(x^2) - xP_1(x ^ 2)
f(x)=P0(x2)+xP1(x2)f(−x)=P0(x2)−xP1(x2)
思考：想快速得到n个点的值，只需要计算得到P0的n/2个的值，同理只需要计算P1的n/2个的值 。f的系数为a0—an，P0的系数为a0，a2…a2m，P1的系数为a0，a2…a2m+1
function f(a, n) :
P0 = f(Aodd, n / 2);
P1 = f(Aeven, n / 2);
merge(P0 + xP1, P0 - xP1);
时间复杂度：O(nlogn)
对于x取值为（1，-1，2，-2，3，-3，4，-4，5，-5）
那么就可以求（1，4，9，16，25）就可以
但是（1，4，9，16，25）就不可以再分了
上面的是对称选择
所以引入下面的复数的办法
开平方。从下面往上看，什么的平方是1，那就是-1和1。那什么的平方是-1呢也就是-1开根号，那就是-i和i。所以第二层有四个值（-i，i，-1，1），那么对这四个值开根号得到
1 / 2 − 1 / 2 i − 1 / 2 + 1 / 2 i − i i − 1 / 2 − 1 / 2 i 1 / 2 + 1 / 2 i − 1 1 \sqrt{1/2}-\sqrt{1/2}i \\ -\sqrt{1/2}+\sqrt{1/2}i\\ -i\\ i\\ -\sqrt{1/2}-\sqrt{1/2}i\\ \sqrt{1/2}+\sqrt{1/2}i\\ -1\\ 1\\ 1/2−1/2i−1/2+1/2i−ii−1/2−1/2i1/2+1/2i−11
复数四则运算：
z1 = a + bi z2 = c + di
加法：z1 + z2 = （a + c） + (b + d)i
减法：z1 - z2 = (a - c) + (b - d)i
乘法：z1 x z2 = (ac - bd) + (ad + bc)i
除法：z1 / z2 = ((ac + bd) + (bc - ad)i ) / c^2 + d^2
a + bi = r1 * (cosA + isinA) c + di = r2 * (cosB + isinB)
(a + bi)(c + di) = r1 * r2 * (cosA + isinA) * (cosB + isinB)
天才想法2-复平面上的单位圆
如果要得到十六个点的值就要对
1
/
2
−
1
/
2
i
−
1
/
2
+
1
/
2
i
−
i
i
−
1
/
2
−
1
/
2
i
1
/
2
+
1
/
2
i
−
1
1
\sqrt{1/2}-\sqrt{1/2}i \\ -\sqrt{1/2}+\sqrt{1/2}i\\ -i\\ i\\ -\sqrt{1/2}-\sqrt{1/2}i\\ \sqrt{1/2}+\sqrt{1/2}i\\ -1\\ 1\\
1/2−1/2i−1/2+1/2i−ii−1/2−1/2i1/2+1/2i−11
八个这样的复杂的点运算，所以这就很复杂，所以引入了下面这种把点放到一个圆上的办法
每个点都是对应a + bi，都有一个b虚部和一个实部a，构成一个坐标（a， b）
以1开始，所以1这个点就是（1，0）
对1开根号得到在圆上1到1的弧的之间那个点-1和与其对称点1
如果对-1开根号就是得到在圆上1到-1的弧上中间点i与其对称点-i
然后对i开根号就得到在圆上1到i的弧的中间那个点
1
/
2
+
1
/
2
i
\sqrt{1/2}+\sqrt{1/2}i\\
1/2+1/2i
和与其对称的点
−
1
/
2
−
1
/
2
i
-\sqrt{1/2}-\sqrt{1/2}i\\
−1/2−1/2i
同理对-i开一次平方得到在圆1到-i的弧之间那个点与其对称的点
1
/
2
−
1
/
2
i
−
1
/
2
+
1
/
2
i
\sqrt{1/2}-\sqrt{1/2}i \\ -\sqrt{1/2}+\sqrt{1/2}i\\
1/2−1/2i−1/2+1/2i
总结：要求某个点N的开根号就是这个点N与1在圆上弧的中间点P和与P的对称点Q
同时发现八个点是圆的八等分点
那么十六个值就是十六等分点
每个点到原点的距离都是1。对于表示标号1这个点的横坐标和纵坐标。那么就可以用角度来表示
对于16个点来说
x = c o s ( 2 Π 16 ) y = s i n ( 2 Π 16 ) x = cos(\frac{2Π}{16})\\ y = sin(\frac{2Π}{16}) x=cos(162Π)y=sin(162Π)
那么就可以得到公式：（n为n等分，k就是第k个点， w就是一个虚数）
w n k = c o s ( 2 Π k n ) + s i n ( 2 Π k n ) i w_n^k = cos(\frac{2Πk}{n}) + sin(\frac{2Πk}{n})i wnk=cos(n2Πk)+sin(n2Πk)i
w相关计算：
( w k n ) 2 = w n 2 k = w n / 2 k (w_k^n)^2 = w^{2k}_n = w^k_{n/2} (wkn)2=wn2k=wn/2k
解释：比如上面那个16等分标号的图，标号4的平方就是标号8。某k位置的点的平方也就是2k位置
−
w
n
k
=
w
n
k
+
n
/
2
=
w
n
/
2
k
-w_n^k = w^{k + n/2}_{n} = w^k_{n / 2}
−wnk=wnk+n/2=wn/2k
解释：比如上面那个16等分标号的图，标号2的对称点就是标号10，标号6的对称点就是标号14，都是相差8。某k位置的点的对称点也就是k + n / 2位置的点
w
n
−
k
=
c
o
s
(
2
Π
k
n
)
−
s
i
n
(
2
Π
k
n
)
i
w^{-k}_n = cos(\frac{2Πk}{n}) - sin(\frac{2Πk}{n})i
wn−k=cos(n2Πk)−sin(n2Πk)i
解释：比如上面那个16等分标号的图，标号15的点其实就是-1点，也就是跟1的纵坐标相反的的·1点。某k位置的点的-k位置的点与k点的纵坐标相反
推断：
根据：
f
(
x
)
=
P
0
(
x
2
)
+
x
P
1
(
x
2
)
f
(
−
x
)
=
P
0
(
x
2
)
−
x
P
1
(
x
2
)
x
=
w
n
k
(
w
k
n
)
2
=
w
n
2
k
=
w
n
/
2
k
−
w
n
k
=
w
n
k
+
n
/
2
=
w
n
/
2
k
得到：
f
(
w
n
k
)
=
P
0
(
w
n
/
2
k
)
+
w
n
k
P
1
(
w
n
/
2
k
)
f
(
w
n
k
+
n
/
2
)
=
P
0
(
w
n
/
2
k
)
−
w
n
k
P
1
(
w
n
/
2
k
)
根据：\\ f(x) = P_0(x^2) + xP_1(x^2)\\ f(-x) = P_0(x^2) - xP_1(x^2)\\ x = w_n^k\\ (w_k^n)^2 = w^{2k}_n = w^k_{n/2}\\ -w_n^k = w^{k + n/2}_{n} = w^k_{n / 2}\\ 得到：\\ f( w_n^k) = P_0(w_{n/2}^k) + w_n^kP_1(w_{n/2}^k)\\ f(w_n^{k+n/2}) = P_0(w_{n/2}^k) - w_n^kP_1(w_{n/2}^k)\\
根据：f(x)=P0(x2)+xP1(x2)f(−x)=P0(x2)−xP1(x2)x=wnk(wkn)2=wn2k=wn/2k−wnk=wnk+n/2=wn/2k得到：f(wnk)=P0(wn/2k)+wnkP1(wn/2k)f(wnk+n/2)=P0(wn/2k)−wnkP1(wn/2k)
FFT算法闪亮登场
f ( w n k ) = P 0 ( w n / 2 k ) + w n k P 1 ( w n / 2 k ) f ( w n k + n / 2 ) = P 0 ( w n / 2 k ) − w n k P 1 ( w n / 2 k ) f( w_n^k) = P_0(w_{n/2}^k) + w_n^kP_1(w_{n/2}^k)\\ f(w_n^{k+n/2}) = P_0(w_{n/2}^k) - w_n^kP_1(w_{n/2}^k) f(wnk)=P0(wn/2k)+wnkP1(wn/2k)f(wnk+n/2)=P0(wn/2k)−wnkP1(wn/2k)
function f(a, n)
P0 = f(Aodd, n/2);
P1 = f(Aeven, n/2);
merger(P0+xP1, P0-P1);
时间复杂度O(nlogn)
代码演示：
系数表示到点值表示
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#include <string>
#include <map>
#include <vector>
#include <set>
#include <stack>
using namespace std;
//复数类
class Complex {
public:
Complex(double r = 0.0, double i = 0.0) : r(r), i(i) {}
Complex &operator*=(const Complex &rhl) {
double a = r, b = i, c = rhl.r, d = rhl.i;
r = a * c - b * d;
i = a * d + b * c;
return *this;
}
Complex operator*(const Complex &rhl) const {
Complex ret(*this);
ret *= rhl;
return ret;
}
Complex &operator+=(const Complex &rhl) {
r += rhl.r;
i += rhl.i;
return *this;
}
Complex operator+(const Complex &rhl) const {
Complex ret(*this);
ret += rhl;
return ret;
}
Complex &operator-=(const Complex &rhl) {
r -= rhl.r;
i -= rhl.i;
return *this;
}
Complex operator-(const Complex &rhl) const {
Complex ret(*this);
ret -= rhl;
return ret;
}
private :
double r, i;
};
struct FastFourierTransform {
#define PI acos(-1)
void transform(vector<Complex> &a, int n) {//一个系数表示集合
if (n == 1) {return ;}
int m = n / 2;
vector<Complex> a1(m), a2(m);
for (int i = 0; i < m; i++) {//把a数组中的系数表示拆分到奇偶数组中
//2 * 0， 2 * 1， 2 * 2...//2 * 0 + 1， 2 * 1 + 1， 2 * 2 + 1...
a1[i] = a[i << 1] , a2[i] = a[i << 1 | 1];
}
//对半分解
transform(a1, m);
transform(a2, m);
Complex wn(1, 0), w1(cos(2 * PI / n), sin(2.0 * PI / n));
for (int i = 0; i < m; i++) {
a[i] = a1[i] + wn * a2[i];
a[i + m] = a1[i] - wn * a2[i];
wn *= w1;//update wn -> w1 = w0 * w1 w2 = w1 * w1
}
return ;
}
#undef PI
} fft;
int main(int argc, char *argv[], char *env[])
{
return 0;
}
实战问题：多项式系数求解
问题描述：给出A（x），B（x）多项式的系数表示，求C（x） = A（x）* B（x）的系数表示
假设：A（x）的系数表示为（1，1），B（x）的系数表示为（1，3）
则：C（x）的系数表示为（1，4，3）
A（x） = x + 1
B（x）= 3x + 1
C（x） = 3x^2 + 4x + 1
注：点值就是点值表示法就是带x的项
所以
A（x）有n项就有n - 1个点值
B（x）有m就有m - 1个点值
C（x）就有n + m - 1项
上面例子，A是两项的，B是两项的，那C一定是三项的。要计算出C的点值出来就要计算出大于或等于三个点值，所以根据上面算法奇偶对半划分，就是要取2的次方，所以就取四个点值
核心思路：系数表示推到点值表示，再由点值表示推到系数表示
1、系数表示法推到点值表示法（得x）
问题描述：给出A（x），B（x）多项式的系数表示，求C（x）=A（x）*B（x）的系数表示
假设：A（x）的系数表示为（1，1），B（x）的系数表示为（1，3）
则：C（x）的系数表示为（1，4，3）
A的4个点值表示：
(
x
0
,
a
0
)
、
(
x
1
,
a
1
)
、
(
x
2
,
a
2
)
、
(
x
3
,
a
3
)
(x_0, a_0)、(x_1, a_1)、(x_2,a_2)、(x_3,a_3)
(x0,a0)、(x1,a1)、(x2,a2)、(x3,a3)
B的4个点值表示：
(
x
0
,
b
0
)
、
(
x
1
,
b
1
)
、
(
x
2
,
b
2
)
、
(
x
3
,
b
3
)
(x_0, b_0)、(x_1, b_1)、(x_2,b_2)、(x_3,b_3)
(x0,b0)、(x1,b1)、(x2,b2)、(x3,b3)
C的4个点值表示：
(
x
0
,
a
0
∗
b
0
)
、
(
x
1
,
a
1
∗
b
1
)
、
(
x
2
,
a
2
∗
b
2
)
、
(
x
3
,
a
3
∗
b
3
)
(x_0, a_0*b_0)、(x_1, a_1*b_1)、(x_2,a_2*b_2)、(x_3,a_3*b_3)
(x0,a0∗b0)、(x1,a1∗b1)、(x2,a2∗b2)、(x3,a3∗b3)
2、由点值表示法返推到系数表示法（求P）
问题描述：给出A（x），B（x)多项式的系数表示，求C（x）=A（x）* B（x）的系数表示
A的4个点值表示：
(
w
4
0
,
a
0
)
、
(
w
4
1
,
a
1
)
、
(
w
4
2
,
a
2
)
、
(
w
4
3
,
a
3
)
(w_4^0,a_0)、(w_4^1,a_1)、(w_4^2,a_2)、(w_4^3,a_3)
(w40,a0)、(w41,a1)、(w42,a2)、(w43,a3)
B的4个点值表示：
(
w
4
0
,
b
0
)
、
(
w
4
1
,
b
1
)
、
(
w
4
2
,
b
2
)
、
(
w
4
3
,
b
3
)
(w_4^0,b_0)、(w_4^1,b_1)、(w_4^2,b_2)、(w_4^3,b_3)
(w40,b0)、(w41,b1)、(w42,b2)、(w43,b3)
C的4个点值表示：
(
w
4
0
,
a
0
∗
b
0
)
、
(
w
4
1
,
a
1
∗
b
1
)
、
(
w
4
2
,
a
2
∗
b
2
)
、
(
w
4
3
,
a
3
∗
b
3
)
(w_4^0,a_0*b_0)、(w_4^1,a_1*b_1)、(w_4^2,a_2*b_2)、(w_4^3,a_3*b_3)
(w40,a0∗b0)、(w41,a1∗b1)、(w42,a2∗b2)、(w43,a3∗b3)
[ P 0 P 1 P 2 P 3 ] = [ 1 w 4 0 ( w 4 0 ) 2 ( w 4 0 ) 3 1 w 4 1 ( w 4 1 ) 2 ( w 4 1 ) 3 1 w 4 2 ( w 4 2 ) 2 ( w 4 2 ) 3 1 w 4 3 ( w 4 3 ) 2 ( w 4 3 ) 3 ] [ C 0 C 1 C 2 C 3 ] \left[\begin{matrix} P_0\\P_1\\P_2\\P_3\\ \end{matrix} \right]= \left[ \begin{matrix} 1 &w_4^0 &(w_4^0)^2 &(w_4^0)^3\\ 1 &w_4^1 &(w_4^1)^2 &(w_4^1)^3\\ 1 &w_4^2 &(w_4^2)^2 &(w_4^2)^3\\ 1 &w_4^3 &(w_4^3)^2 &(w_4^3)^3\\ \end{matrix} \right] \left[\begin{matrix} C_0\\ C_1\\ C_2\\ C_3\\ \end{matrix} \right] P0P1P2P3 = 1111w40w41w42w43(w40)2(w41)2(w42)2(w43)2(w40)3(w41)3(w42)3(w43)3 C0C1C2C3
P = W C W − 1 P = W − 1 W C = E C = C P = WC\\ W^{-1}P = W^{-1}WC = EC = C P=WCW−1P=W−1WC=EC=C
求得：
[
1
/
4
w
4
0
/
4
(
w
4
0
)
2
/
4
(
w
4
0
)
3
/
4
1
/
4
w
4
−
1
/
4
(
w
4
−
1
)
2
/
4
(
w
4
−
1
)
3
/
4
1
/
4
w
4
−
2
/
4
(
w
4
−
2
)
2
/
4
(
w
4
−
2
)
3
/
4
1
/
4
w
4
−
3
/
4
(
w
4
−
3
)
2
/
4
(
w
4
−
3
)
3
/
4
]
[
P
0
P
1
P
2
P
3
]
=
[
C
0
C
1
C
2
C
3
]
\left[ \begin{matrix} 1/4 &w_4^0/4 &(w_4^0)^2/4 &(w_4^0)^3/4\\ 1/4 &w_4^{-1}/4 &(w_4^{-1})^2/4 &(w_4^{-1})^3/4\\ 1/4 &w_4^{-2}/4 &(w_4^{-2})^2/4 &(w_4^{-2})^3/4\\ 1/4 &w_4^{-3}/4 &(w_4^{-3})^2/4 &(w_4^{-3})^3/4\\ \end{matrix} \right] \left[\begin{matrix} P_0\\P_1\\P_2\\P_3\\ \end{matrix} \right]= \left[\begin{matrix} C_0\\ C_1\\ C_2\\ C_3\\ \end{matrix} \right]
1/41/41/41/4w40/4w4−1/4w4−2/4w4−3/4(w40)2/4(w4−1)2/4(w4−2)2/4(w4−3)2/4(w40)3/4(w4−1)3/4(w4−2)3/4(w4−3)3/4
P0P1P2P3
=
C0C1C2C3
1 4 [ 1 w 4 0 ( w 4 0 ) 2 ( w 4 0 ) 3 1 w 4 − 1 ( w 4 − 1 ) 2 ( w 4 − 1 ) 3 1 w 4 − 2 ( w 4 − 2 ) 2 ( w 4 − 2 ) 3 1 w 4 − 3 ( w 4 − 3 ) 2 ( w 4 − 3 ) 3 ] [ P 0 P 1 P 2 P 3 ] = [ C 0 C 1 C 2 C 3 ] \frac{1}{4}\left[ \begin{matrix} 1 &w_4^0 &(w_4^0)^2 &(w_4^0)^3\\ 1 &w_4^{-1} &(w_4^{-1})^2 &(w_4^{-1})^3\\ 1 &w_4^{-2} &(w_4^{-2})^2 &(w_4^{-2})^3\\ 1 &w_4^{-3} &(w_4^{-3})^2 &(w_4^{-3})^3\\ \end{matrix} \right] \left[\begin{matrix} P_0\\P_1\\P_2\\P_3\\ \end{matrix} \right]= \left[\begin{matrix} C_0\\ C_1\\ C_2\\ C_3\\ \end{matrix} \right] 41 1111w40w4−1w4−2w4−3(w40)2(w4−1)2(w4−2)2(w4−3)2(w40)3(w4−1)3(w4−2)3(w4−3)3 P0P1P2P3 = C0C1C2C3
补充逆矩阵及其求解办法：
若AB = BA = E，则称A为可逆矩阵
且存在性质：
1、若存在AB = BA = E, AC = CA = E。就存在B = BE = B(AC) = (BA)C = EC = C
2、AA* = |A|E A为方阵。|A|为行列式。A*为转置矩阵，也就是代数余子式集合2
求解逆矩阵：
1、AA* = |A|E
2、初等变换
（1）构造nx2n维矩阵（A ：E）
（2）对（A：E）施行一系列初等行变换，直至将其左边子矩阵A化为单位矩阵E，此时右边子矩阵即为$ A_-1 $
（3）待定系数法。
代码演示：
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <queue>
#include <stack>
#include <algorithm>
#include <string>
#include <map>
#include <set>
#include <vector>
using namespace std;
class Complex {
public :
Complex(double r = 0.0, double i = 0.0) : r(r), i(i) {}
Complex &operator*=(const Complex &rhl) {
double a = r, b = i, c = rhl.r, d = rhl.i;
r = a * c - b * d;
i = a * d + b * c;
return *this;
}
Complex operator*(const Complex &rhl) const {
Complex ret(*this);
ret *= rhl;
return ret;
}
Complex &operator+=(const Complex &rhl) {
r += rhl.r;
i += rhl.i;
return *this;
}
Complex operator+(const Complex &rhl) const {
Complex ret(*this);
ret += rhl;
return ret;
}
Complex &operator-=(const Complex &rhl) {
r -= rhl.r;
i -= rhl.i;
return *this;
}
Complex operator-(const Complex &rhl) const {
Complex ret(*this);
ret -= rhl;
return ret;
}
Complex &operator/=(const double x) {
r /= x;
i /= x;
return *this;
}
Complex operator/(const double x) const {
Complex ret(*this);
ret /= x;
return ret;
}
double r, i;
};
ostream &operator<<(ostream &out, const Complex &a) {
out << a.r << "+" << a.i << "i";
return out;
}
struct FastFourierTransform {
#define PI acos(-1)
void __transform(vector<Complex> &a, int n, int type = 1) {//type代表sin前面得符号
if (n == 1) { return ; }
int m = n / 2;
vector<Complex> a1(m), a2(m);
for (int i = 0; i < m; i++) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
__transform(a1, m, type);
__transform(a2, m, type);
Complex wn(1, 0), w1(cos(2.0 * PI / n), type * sin(2.0 * PI / n));
for (int i = 0; i < m; i++) {
a[i] = a1[i] + wn * a2[i];
a[i + m] = a1[i] - wn * a2[i];
wn *= w1;
}
return ;
}
//得到点值表示
vector<Complex> dft(vector<Complex> &a, int n) {
vector<Complex> temp(n);
for (int i = 0, m = a.size(); i < m; i++) temp[i] = a[i];//保存a
__transform(temp, n);
return temp;
}
//得到系数表示
vector<Complex> idft(vector<Complex> &a, int n) {
vector<Complex> temp(n);
for (int i = 0, m = a.size(); i < m; i++) temp[i] = a[i];
__transform(temp, n, -1);
for (int i = 0; i < n; i++) temp[i] /= n;
return temp;
}
#undef PI
} fft;
int main() {
int n, m, N = 1;
cin >> n >> m;//A多项式有n项，B多项式有m项
while (N < n + m - 1) N *= 2;//求得C的最小二倍项数
vector<Complex> a(n), b(m);//存的是AB多项式的系数表示法
for (int i = 0; i < n; i++) cin >> a[i].r;
for (int i = 0; i < m; i++) cin >> b[i].r;
vector<Complex> a_val = fft.dft(a, N);//求A多项式的点值表示
vector<Complex> b_val = fft.dft(b, N);//求B多项式的点值表示
vector<Complex> c_val(N);//
for (int i = 0; i < N; i++) c_val[i] = a_val[i] * b_val[i];//得到C多项式的点值表示法
for (int i = 0; i < N; i++) cout << c_val[i] << "\t";
cout << endl;
vector<Complex> c = fft.idft(c_val, N);//得到C多项式的系数表示法
for (int i = 0; i < n + m - 1; i++) cout << c[i].r << " ";
cout << endl;
return 0;
}
声音处理
认识声音
认识声音：
首先要感受1Hz到5Hz正弦波的声音（声音文件私聊）
听完会发现从1Hz到5Hz越来越快
1Hz代表正弦波在1秒内走过1个周期
2Hz代表正弦波在1秒内走过2个周期
3Hz代表正弦波在1秒内走过3个周期
4Hz代表正弦波在1秒内走过4个周期
5Hz代表正弦波在1秒内走过5个周期
声音在计算机中存储：
采样频率：5Hz -> 1秒采样5个点
声音采样点图：（此图由py文件绘制生成）
这个图看起来是连续的，但是不是连续的，只是因为采样频率大就看起来就连续
制作声音文件：
#!/usr/bin/env python
# coding=utf-8
import wave
import numpy as np
frameRate = 16100
time = 10
volumn = 40
pi = np.pi
t = np.arange(0, time, 1.0 / frameRate) #从0开始生成time个值，每两个值之间相差1.0 / frameRate
def gen_wave_data(fileName, readlf):
wave_data = volumn * np.sin(2.0 * pi * readlf * t);
wave_data = wave_data.astype(np.int8)
f = wave.open(fileName, "wb")
f.setnchannels(1) #单声道，双声道使声音变立体
f.setsampwidth(1) #设置采样的位宽，也就是占用内存，比如无损音乐的位宽至少2字节，也就是0~65535
f.setframerate(frameRate) #告诉计算机一秒多少个数据
f.writeframes(wave_data.tostring())
f.close();
print("write data to : " + fileName)
gen_wave_data("1Hz.wav", 1)
gen_wave_data("2Hz.wav", 2)
gen_wave_data("3Hz.wav", 3)
gen_wave_data("4Hz.wav", 4)
gen_wave_data("5Hz.wav", 5)
绘制声音文件：
#!/usr/bin/env python
# coding=utf-8
import wave
import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as pl
import numpy as np
def draw_wava_data(fileName, plotNum):
f = wave.open(fileName, "rb")
params = f.getparams()
str_data = f.readframes(params[3])
f.close()
wave_data = np.fromstring(str_data, dtype = np.int8)
t = np.arange(0, len(wave_data) * 1.0 / params[2], 1.0 / params[2])
pl.subplot(plotNum)
pl.title(fileName)
pl.plot(t, wave_data)
pl.xlabel("time(seconds)")
draw_wava_data("1Hz.wav", 321)
draw_wava_data("2Hz.wav", 322)
draw_wava_data("3Hz.wav", 323)
draw_wava_data("4Hz.wav", 324)
draw_wava_data("5Hz.wav", 325)
pl.show()
DFT算法简介(离散傅里叶变换)
作用：周期信号分解（傅里叶老头子认为所有的信号都是呈周期性，所以也可以称为信号分解）
比如一个复杂信号是由1Hz和2Hz组成的，DFT算法就是把这个复杂信号分解为1Hz和2Hz
公式：
X
(
m
)
=
∑
n
=
0
N
−
1
x
(
n
)
⋅
e
−
j
2
Π
m
n
/
N
X(m) = \sum_{n = 0}^{N - 1}{x(n) · e ^ {-j2Πmn/N}}
X(m)=n=0∑N−1x(n)⋅e−j2Πmn/N
公式解析：X(m)代表第m种频率的分析结果，x（n）中的n代表第n时间点的采样信号数据，对于整个e中的m代表m种频率
补充基础知识： e j 2 Π k / n = w n k = c o s ( 2 Π k n ) + s i n ( 2 Π k n ) i e − j 2 Π k / n = w n − k = c o s ( 2 Π k n ) − s i n ( 2 Π k n ) i 补充基础知识：\\ e^{j2Πk/n} = w^k_n = cos(\frac{2Πk}{n}) + sin(\frac{2Πk}{n})i\\ e^{-j2Πk/n} = w^{-k}_n = cos(\frac{2Πk}{n}) - sin(\frac{2Πk}{n})i 补充基础知识：ej2Πk/n=wnk=cos(n2Πk)+sin(n2Πk)ie−j2Πk/n=wn−k=cos(n2Πk)−sin(n2Πk)i
得到： X ( m ) = ∑ N = 0 N − 1 x ( n ) ⋅ ( w N − m ) n 得到：\\ X(m) = \sum_{N=0}^{N - 1}x(n)·(w^{-m}_{N})^n 得到：X(m)=N=0∑N−1x(n)⋅(wN−m)n
[ 1 w 4 0 ( w 4 0 ) 2 ( w 4 0 ) 3 1 w 4 − 1 ( w 4 − 1 ) 2 ( w 4 − 1 ) 3 1 w 4 − 2 ( w 4 − 2 ) 2 ( w 4 − 2 ) 3 1 w 4 − 3 ( w 4 − 3 ) 2 ( w 4 − 3 ) 3 ] [ x ( 0 ) x ( 1 ) x ( 2 ) x ( 3 ) ] = [ X ( 0 ) X ( 1 ) X ( 2 ) X ( 3 ) ] \left[\begin{matrix} 1 &w_4^0 &(w_4^0)^2 &(w_4^0)^3\\ 1 &w_4^{-1} &(w_4^{-1})^2 &(w_4^{-1})^3\\ 1 &w_4^{-2} &(w_4^{-2})^2 &(w_4^{-2})^3\\ 1 &w_4^{-3} &(w_4^{-3})^2 &(w_4^{-3})^3\\ \end{matrix} \right] \left[\begin{matrix} x(0)\\ x(1)\\ x(2)\\ x(3)\\ \end{matrix} \right]= \left[\begin{matrix} X(0)\\ X(1)\\ X(2)\\ X(3) \end{matrix}\right] 1111w40w4−1w4−2w4−3(w40)2(w4−1)2(w4−2)2(w4−3)2(w40)3(w4−1)3(w4−2)3(w4−3)3 x(0)x(1)x(2)x(3) = X(0)X(1)X(2)X(3)
FFT与DFT对比：
FFT代表一种数据处理的办法，就是一种算法
DFT代表具体一类问题，用来分解信号的，把时域数据（根据时间轴顺序给出的数据）转换到频域数据，比如上面矩阵的x(n) 就是时域，X(m)就是频域
实战1：声音还原
思路：解析声音数据->DFT算法->分析频域数据->还原声音
1、解析声音数据：（会产生一个一串数字文件）
#!/usr/bin/env python
# coding=utf-8
import wave
import numoty as np
f = wave.open("xxx.wave", "rb");//xxx.wave为一个复杂声音文件
params = f.getparams()
str_data = f.readframes(params[3])
f.close()
wave_data = np.fromstring(str_data, dtype = np.int8)
output_fileName = "input.data"
fout = open(output_fileName, "w")
fout.write(str(params[3] + "\n"))
for x in wave_data:
fout.write(str(x) + "\n")
fout.close()
print("write wave to : " + output_fileName)
2、DFT算法：
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <queue>
#include <stack>
#include <algorithm>
#include <string>
#include <map>
#include <set>
#include <vector>
using namespace std;
class Complex {
public :
Complex(double r = 0.0, double i = 0.0) : r(r), i(i) {}
Complex &operator*=(const Complex &rhl) {
double a = r, b = i, c = rhl.r, d = rhl.i;
r = a * c - b * d;
i = a * d + b * c;
return *this;
}
double m() {return sqrt(r * r + i * i);}
Complex operator*(const Complex &rhl) const {
Complex ret(*this);
ret *= rhl;
return ret;
}
Complex &operator+=(const Complex &rhl) {
r += rhl.r;
i += rhl.i;
return *this;
}
Complex operator+(const Complex &rhl) const {
Complex ret(*this);
ret += rhl;
return ret;
}
Complex &operator-=(const Complex &rhl) {
r -= rhl.r;
i -= rhl.i;
return *this;
}
Complex operator-(const Complex &rhl) const {
Complex ret(*this);
ret -= rhl;
return ret;
}
Complex &operator/=(const double x) {
r /= x;
i /= x;
return *this;
}
Complex operator/(const double x) const {
Complex ret(*this);
ret /= x;
return ret;
}
double r, i;
};
ostream &operator<<(ostream &out, const Complex &a) {
out << a.r << "+" << a.i << "i";
return out;
}
struct FastFourierTransform {
#define PI acos(-1)
void __transform(vector<Complex> &a, int n, int type = 1) {//type代表sin前面得符号
if (n == 1) { return ; }
int m = n / 2;
vector<Complex> a1(m), a2(m);
for (int i = 0; i < m; i++) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
__transform(a1, m, type);
__transform(a2, m, type);
Complex wn(1, 0), w1(cos(2.0 * PI / n), type * sin(2.0 * PI / n));
for (int i = 0; i < m; i++) {
a[i] = a1[i] + wn * a2[i];
a[i + m] = a1[i] - wn * a2[i];
wn *= w1;
}
return ;
}
//得到点值表示
vector<Complex> dft(vector<Complex> &a, int n) {
vector<Complex> temp(n);
for (int i = 0, m = a.size(); i < m; i++) temp[i] = a[i];//保存a
__transform(temp, n, -1);//跟fft基础上换为负
return temp;
}
//得到系数表示
vector<Complex> idft(vector<Complex> &a, int n) {
vector<Complex> temp(n);
for (int i = 0, m = a.size(); i < m; i++) temp[i] = a[i];
__transform(temp, n);//在fft基础上换为正
for (int i = 0; i < n; i++) temp[i] /= n;
return temp;
}
#undef PI
} fft;
int main() {
int n, N = 1;
cin >> n;
vector<Complex> x(n);
for (int i = 0; i < n; i++) cin >> x[i].r;
while (N < n) N <<= 1;
vector<Complex> X = fft.dft(x, N);
for (int i = 0; i < N; i++) {
cout << X[i].r << " " << X[i].i << " " << X[i].m() / N << endl;
}
return 0;
}
3、分析频域数据
对上面离散傅里叶变换得到的结果转换的图
对于上面的图，根据DFT算法性质是对称的，所以只有一半的图是有用的。只要看柱的高度就可以，柱的高度代表声音的大小，底是频率
4、还原声音
#!/usr/bin/env python
# coding=utf-8
import wave
import numpy as np
frameRate = 16100
time = 10
volumn = 40
pi = np.pi
t = np.arange(0, time, 1.0 / frameRate) #从0开始生成time个值，每两个值之间相差1.0 / frameRate
fileName = "guess.wav"
#wave_data = volumn * np.sin(2.0 * pi * readlf * t);#readilf为真实频域，volumn为声音大小，这两个就是通过DFT得到的
wave_data = 30 * 14 * np.sin(2.0 * pi * 50 * t) + np.sin(2.0 * pi * 2000 * t)
wave_data = wave_data.astype(np.int8)
f = wave.open(fileName, "wb")
f.setnchannels(1) #单声道，双声道使声音变立体
f.setsampwidth(1) #设置采样的位宽，也就是占用内存，比如无损音乐的位宽至少2字节，也就是0~65535
f.setframerate(frameRate) #告诉计算机一秒多少个数据
f.writeframes(wave_data.tostring())
f.close();
print("write data to : " + fileName)
实战2：声音减噪
一段复杂声音结果DFT转换得到的图像:
思考：什么频率的信号是噪音呢？声音比较小频率比较高的声音，因为太高频率的频率的声音对于我们人来说是听不到的，所以没有什么意义。所以降噪就是把高频信号全部过滤掉。
因为DFT得到的结果是对称的，所以上面图越接近中间点位置就越是高频信号
拓展：FFT算法可以用在信号处理等场景，比如PS抠图
实战未完待续
500
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
