[数据结构与算法]傅里叶变换-信号处理

傅里叶变换

系数表示法与点值表示法

$$f(x)=a_0?x^0+a_1?x^1+a_2?x^2+a_3?x^3+...+a_n?x^n$$

系数表示法：
$$(a_0,a_1,a_2,a_3,a_4,...,a_n)$$

$$例如(1,2,1,4)=>f(x)=1+2?x^1+1?x^2+4?x^3$$

点值表示法：n + 1个点坐标确定一条曲线
$$(x_0,y_0)、(x_1,y_1)、(x_2,y_2)、(x_3,y_3).....、(x_n,y_n)$$

FFT算法的作用

$$\begin{gathered} 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) \\ \sum _0^{n}f(x) \end{gathered}$$

等价于：

$$?\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ .\\ .\\ .\\ y_n\end{array}?=?\begin{array}{cccccc}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\end{array}?\times ?\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ .\\ .\\ .\\ a_n\end{array}?$$

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$$

思考：在一个多项式里面是否存在对称性，如果存在就只要求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)$$

第二、多项式构造
$$\begin{gathered} 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) \end{gathered}$$
第三、抽象成算法
$$\begin{gathered} f(x)=P_0(x^2)+x{P}_1(x^2) \\ f(?x)=P_0(x^2)?x{P}_1(x^2) \end{gathered}$$

思考：想快速得到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），那么对这四个值开根号得到

$$\begin{gathered} \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 \end{gathered}$$

复数四则运算：

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-复平面上的单位圆

如果要得到十六个点的值就要对
$$\begin{gathered} \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 \end{gathered}$$
八个这样的复杂的点运算，所以这就很复杂，所以引入了下面这种把点放到一个圆上的办法

每个点都是对应a + bi，都有一个b虚部和一个实部a，构成一个坐标（a， b）

以1开始，所以1这个点就是（1，0）

对1开根号得到在圆上1到1的弧的之间那个点-1和与其对称点1

如果对-1开根号就是得到在圆上1到-1的弧上中间点i与其对称点-i

然后对i开根号就得到在圆上1到i的弧的中间那个点
$$\sqrt{1/2}+\sqrt{1/2}i$$
和与其对称的点
$$?\sqrt{1/2}?\sqrt{1/2}i$$

同理对-i开一次平方得到在圆1到-i的弧之间那个点与其对称的点
$$\begin{gathered} \sqrt{1/2}?\sqrt{1/2}i \\ ?\sqrt{1/2}+\sqrt{1/2}i \end{gathered}$$

总结：要求某个点N的开根号就是这个点N与1在圆上弧的中间点P和与P的对称点Q

同时发现八个点是圆的八等分点
那么十六个值就是十六等分点

每个点到原点的距离都是1。对于表示标号1这个点的横坐标和纵坐标。那么就可以用角度来表示

对于16个点来说

$$\begin{gathered} x=cos(\frac{2\Pi }{16}) \\ y=sin(\frac{2\Pi }{16}) \end{gathered}$$

那么就可以得到公式：（n为n等分，k就是第k个点， w就是一个虚数）

$$w_n^k=cos(\frac{2\Pi k}{n})+sin(\frac{2\Pi k}{n})i$$

w相关计算：

$$(w_k^n{)}^2=w_n^{2k}=w_{n/2}^k$$

解释：比如上面那个16等分标号的图，标号4的平方就是标号8。某k位置的点的平方也就是2k位置

$$?w_n^k=w_n^{k+n/2}=w_{n/2}^k$$
解释：比如上面那个16等分标号的图，标号2的对称点就是标号10，标号6的对称点就是标号14，都是相差8。某k位置的点的对称点也就是k + n / 2位置的点
$$w_n^{?k}=cos(\frac{2\Pi k}{n})?sin(\frac{2\Pi k}{n})i$$
解释：比如上面那个16等分标号的图，标号15的点其实就是-1点，也就是跟1的纵坐标相反的的·1点。某k位置的点的-k位置的点与k点的纵坐标相反

推断：
$$\begin{gathered} 根据： \\ 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^{2k}=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) \end{gathered}$$

FFT算法闪亮登场

$$\begin{gathered} 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) \end{gathered}$$

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)$$
B的4个点值表示：
$$(x_0,b_0)、(x_1,b_1)、(x_2,b_2)、(x_3,b_3)$$
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)$$

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)$$
B的4个点值表示：
$$(w_4^0,b_0)、(w_4^1,b_1)、(w_4^2,b_2)、(w_4^3,b_3)$$
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)$$

$$?\begin{array}{c}{P}_0\\ P_1\\ P_2\\ P_3\end{array}?=?\begin{array}{cccc}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{array}??\begin{array}{c}{C}_0\\ C_1\\ C_2\\ C_3\end{array}?$$

$$\begin{gathered} P=WC \\ {W}^{?1}P=W^{?1}WC=EC=C \end{gathered}$$

求得：
$$?\begin{array}{cccc}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{array}??\begin{array}{c}{P}_0\\ P_1\\ P_2\\ P_3\end{array}?=?\begin{array}{c}{C}_0\\ C_1\\ C_2\\ C_3\end{array}?$$

$$\frac{1}{4}?\begin{array}{cccc}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{array}??\begin{array}{c}{P}_0\\ P_1\\ P_2\\ P_3\end{array}?=?\begin{array}{c}{C}_0\\ C_1\\ C_2\\ C_3\end{array}?$$

补充逆矩阵及其求解办法：

若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)=\sum _{n=0}^{N?1}x(n)?e^{?j2\Pi mn/N}$$

公式解析：X(m)代表第m种频率的分析结果，x（n）中的n代表第n时间点的采样信号数据，对于整个e中的m代表m种频率

$$\begin{gathered} 补充基础知识： \\ {e}^{j2\Pi k/n}=w_n^k=cos(\frac{2\Pi k}{n})+sin(\frac{2\Pi k}{n})i \\ {e}^{?j2\Pi k/n}=w_n^{?k}=cos(\frac{2\Pi k}{n})?sin(\frac{2\Pi k}{n})i \end{gathered}$$

$$\begin{gathered} 得到： \\ X(m)=\sum _{N=0}^{N?1}x(n)?(w_N^{?m}{)}^n \end{gathered}$$

$$?\begin{array}{cccc}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{array}??\begin{array}{c}x(0)\\ x(1)\\ x(2)\\ x(3)\end{array}?=?\begin{array}{c}X(0)\\ X(1)\\ X(2)\\ X(3)\end{array}?$$

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抠图

实战未完待续

开发者涨薪指南

48位大咖的思考法则、工作方式、逻辑体系