Fourier transform signal processing

The Fourier transform

Coefficient representation and point value representation

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

Coefficient representation ：
$$(a_0,a_1,a_2,a_3,a_4,...,a_n)$$

$$exampleSuch\; as(1,2,1,4)=>f(x)=1+2\ast x^1+1\ast x^2+4\ast x^3$$

Point value representation ：n + 1 Point coordinates determine a curve
$$(x_0,y_0)、(x_1,y_1)、(x_2,y_2)、(x_3,y_3).....、(x_n,y_n)$$

FFT The function of algorithm

$$\begin{gathered} f(x)=a_0+a_1\ast x_n^1+a_2\ast x_n^2+a_3\ast x_n^3+\cdot \cdot \cdot +a_m\ast x_n^m(n=1,2,3....m) \\ \sum _0^{n}f(x) \end{gathered}$$

Equivalent to ：

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

FFT The algorithm is to calculate this formula quickly , from n A coefficient is obtained n A value constitutes a point value that represents

Genius idea 1- Quick value

reflection ： How to get a value of six points in a binomial

explain ： Calculate three times to get three points , Six points are obtained by taking symmetric values

$$f(x)=a_0\ast x^0+a_1\ast x^1+a_2\ast x^2+a_3\ast x^3+a_4\ast x^4+a_5\ast x^5$$

reflection ： Whether there is symmetry in a polynomial , If it exists, it only requires n / 2 A little bit

explain ：

First of all 、 Parity partition

$$f(x)=(a_0\ast x^0+a_2\ast x^2+a_4\ast x^4)+(a_1\ast x^1+a_3\ast x^3+a_5\ast x^5)$$

second 、 Polynomial construction
$$\begin{gathered} P_0(y=x^2)=(a_0+a_2\ast y+a_4\ast y) \\ {P}_1(y=x^2)=(a_1+a_3\ast y+a_5\ast y^2) \\ f(x)=P_0(x^2)+x{P}_1(x^2) \\ have\; to：f(-x)=P_0(x^2)-x{P}_1(x^2) \end{gathered}$$
Third 、 Abstract it into an algorithm
$$\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}$$

reflection ： Want to get... Quickly n The value of a point , Just calculate P0 Of n/2 The value of , Similarly, we only need to calculate P1 Of n/2 The value of .f The coefficient of is a0—an,P0 The coefficient of is a0,a2…a2m,P1 The coefficient of is a0,a2…a2m+1

function f(a, n) :
	P0 = f(Aodd, n / 2);
	P1 = f(Aeven, n / 2);
	merge(P0 + xP1, P0 - xP1);

Time complexity ：O(nlogn)

about x The value is （1,-1,2,-2,3,-3,4,-4,5,-5）

Then you can ask （1,4,9,16,25） Can

however （1,4,9,16,25） You can't divide any more

Above is a symmetric selection

So introduce the following plural method

Square root . Look up from below , The square of what is 1, That's it -1 and 1. So the square of something is -1 That is -1 Square root , That's it -i and i. So the second layer has four values （-i,i,-1,1）, Then we can take the root of these four values and get

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

Four operations of complex number ：

z1 = a + bi z2 = c + di

Add ：z1 + z2 = （a + c） + (b + d)i

Subtraction ：z1 - z2 = (a - c) + (b - d)i

Multiplication ：z1 x z2 = (ac - bd) + (ad + bc)i

division ：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)

Genius idea 2- Unit circle on complex plane

If you want to get the value of sixteen points, you have to be right
$$\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}$$
Eight such complex point operations , So it's complicated , So the following method of putting points on a circle is introduced

Each point corresponds to a + bi, There is one. b The imaginary part and a real part a, Form a coordinate （a, b）

With 1 Start , therefore 1 This point is （1,0）

Yes 1 The open root sign gets on the circle 1 To 1 The point between the arcs of -1 And its point of symmetry 1

If the -1 To open the root sign is to get on the circle 1 To -1 The middle point on the arc of i Its point of symmetry -i

Then on i Open the root sign to get on the circle 1 To i The point in the middle of the arc of
$$\sqrt{1/2}+\sqrt{1/2}i$$
And a point symmetrical to it
$$-\sqrt{1/2}-\sqrt{1/2}i$$

Same reason -i Open the square once to get in the circle 1 To -i The point between the arcs of and its symmetrical point
$$\begin{gathered} \sqrt{1/2}-\sqrt{1/2}i \\ -\sqrt{1/2}+\sqrt{1/2}i \end{gathered}$$

summary ： Ask for a point N The open radical of is this point N And 1 At the middle point of the arc on the circle P And P The point of symmetry of Q

At the same time, it is found that the eight points are the octant of the circle
So the sixteen values are the sixteenth equinox

The distance from each point to the origin is 1. For the representation label 1 The abscissa and ordinate of this point . Then it can be expressed in terms of angle

about 16 For one thing

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

Then you can get the formula ：（n by n Equal division ,k That is the first. k A little bit , w Is an imaginary number ）

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

w Relevant calculations ：

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

explain ： Like the one above 16 A graph of equal division marks , label 4 The square of is the label 8. some k The square of the point of position is 2k Location

$$-w_n^k=w_n^{k+n/2}=w_{n/2}^k$$
explain ： Like the one above 16 A graph of equal division marks , label 2 The symmetry point of is the label 10, label 6 The symmetry point of is the label 14, It's all a difference 8. some k The point of symmetry of the point of position is k + n / 2 The point of location
$$w_n^{-k}=cos(\frac{2\Pi k}{n})-sin(\frac{2\Pi k}{n})i$$
explain ： Like the one above 16 A graph of equal division marks , label 15 The point is -1 spot , That is to say, with 1 The vertical coordinates of the are opposite ·1 spot . some k The position of the point of -k The point of position and k The ordinates of the points are opposite

infer ：
$$\begin{gathered} rootAccording\; to\; the： \\ 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 \\ have\; toTo： \\ 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 Algorithms make their debut

$$\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);

Time complexity O(nlogn)

Code demonstration ：

The coefficient represents to the point value represents

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#include <string>
#include <map>
#include <vector>
#include <set>
#include <stack>
using namespace std;

// Plural 
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) {
    // A coefficient represents a set 
        if (n == 1) {
    return ;}
        int  m = n / 2;
        vector<Complex> a1(m), a2(m);
        for (int i = 0; i < m; i++) {
    // hold a The coefficients in the array represent splitting into odd and even arrays 
            //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];
            
        }
        // Split in half 
        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;
}

Practical problems ： Polynomial coefficient solution

Problem description ： give A（x）,B（x） The coefficient representation of a polynomial , seek C（x） = A（x）* B（x） Coefficient expression of

hypothesis ：A（x） The coefficient of is expressed as （1,1）,B（x） The coefficient of is expressed as （1,3）

be ：C（x） The coefficient of is expressed as （1,4,3）

A（x） = x + 1

B（x）= 3x + 1

C（x） = 3x^2 + 4x + 1

notes ： A point value is a point value, and its representation is a band x The item

therefore

A（x） Yes n There are items n - 1 A point value

B（x） Yes m There is m - 1 A point value

C（x） There is n + m - 1 term

The above example ,A It's a two pronged ,B It's a two pronged , that C It must be trinomial . To work out C The point value of is greater than or equal to three points , So according to the above algorithm, odd even half Division , Is to take 2 To the power of , So take four values

The core idea ： The coefficient represents pushing to the point value represents , Then from the point value representation to the coefficient representation

1、 The coefficient representation is pushed to the point value representation （ have to x）

Problem description ： give A（x）,B（x） The coefficient representation of a polynomial , seek C（x）=A（x）*B（x） Coefficient expression of

hypothesis ：A（x） The coefficient of is expressed as （1,1）,B（x） The coefficient of is expressed as （1,3）

be ：C（x） The coefficient of is expressed as （1,4,3）

A Of 4 A point value indicates ：
$$(x_0,a_0)、(x_1,a_1)、(x_2,a_2)、(x_3,a_3)$$
B Of 4 A point value indicates ：
$$(x_0,b_0)、(x_1,b_1)、(x_2,b_2)、(x_3,b_3)$$
C Of 4 A point value indicates ：
$$(x_0,a_0\ast b_0)、(x_1,a_1\ast b_1)、(x_2,a_2\ast b_2)、(x_3,a_3\ast b_3)$$

2、 From point value representation to coefficient representation （ seek P）

Problem description ： give A（x）,B（x) The coefficient representation of a polynomial , seek C（x）=A（x）* B（x） Coefficient expression of

A Of 4 A point value indicates ：
$$(w_4^0,a_0)、(w_4^1,a_1)、(w_4^2,a_2)、(w_4^3,a_3)$$
B Of 4 A point value indicates ：
$$(w_4^0,b_0)、(w_4^1,b_1)、(w_4^2,b_2)、(w_4^3,b_3)$$
C Of 4 A point value indicates ：
$$(w_4^0,a_0\ast b_0)、(w_4^1,a_1\ast b_1)、(w_4^2,a_2\ast b_2)、(w_4^3,a_3\ast b_3)$$

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

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

Get ：
$$\left[\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}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ P_2\\ P_3\end{array}\right]=\left[\begin{array}{c}{C}_0\\ C_1\\ C_2\\ C_3\end{array}\right]$$

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

Supplementary inverse matrix and its solution ：

if AB = BA = E, said A Is an invertible matrix

And the existence of nature ：

1、 If exist AB = BA = E, AC = CA = E. There is a B = BE = B(AC) = (BA)C = EC = C

2、AA* = |A|E A For the square .|A| Is a determinant .A* For transpose matrix , That is, the set of algebraic cofactors 2

Solving the inverse matrix ：

1、AA* = |A|E

2、 Elementary transformation

（1） structure nx2n D matrix （A ：E）

（2） Yes （A：E） Perform a series of elementary line transformations , Until its left submatrix A Into an identity matrix E, At this point, the right submatrix is $ A_-1 $

（3） Undetermined coefficient method .

Code demonstration ：

#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 representative sin The preceding symbol 
        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 ;
    }
    // Get a point value to represent 
    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];// preservation a
        __transform(temp, n);
        return temp;
    }
    // Get a coefficient to represent 
    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 Polynomials have n term ,B Polynomials have m term 
    while (N < n + m - 1) N *= 2;// Get C The minimum number of binomials of 
    vector<Complex> a(n), b(m);// Deposit is AB Coefficient representation of polynomials  
    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);// seek A Point valued representation of polynomials 
    vector<Complex> b_val = fft.dft(b, N);// seek B Point valued representation of polynomials 
    vector<Complex> c_val(N);//
    for (int i = 0; i < N; i++) c_val[i] = a_val[i] * b_val[i];// obtain C Point valued representation of polynomials 
    for (int i = 0; i < N; i++) cout << c_val[i] << "\t";
    cout << endl;
    vector<Complex> c = fft.idft(c_val, N);// obtain C Coefficient representation of polynomials 
    for (int i = 0; i < n + m - 1; i++) cout << c[i].r << " ";
    cout << endl;
    return 0;
}

Sound processing

Recognize sound

Recognize sound ：

First of all, feel 1Hz To 5Hz The sound of a sine wave （ Voice file private chat ）

After listening, you will find that 1Hz To 5Hz Faster and faster

1Hz Represents the sine wave at 1 In seconds 1 A cycle

2Hz Represents the sine wave at 1 In seconds 2 A cycle

3Hz Represents the sine wave at 1 In seconds 3 A cycle

4Hz Represents the sine wave at 1 In seconds 4 A cycle

5Hz Represents the sine wave at 1 In seconds 5 A cycle

Sound is stored in the computer ：

sampling frequency ：5Hz -> 1 Second sampling 5 A little bit

Sound sampling point diagram ：（ This drawing is provided by py File drawing generation ）

This graph seems to be continuous , But not continuous , Just because the sampling frequency is large, it looks continuous

Making sound files ：

#!/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) # from 0 To start generating time It's worth , The difference between every two values is 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) # Mono , Two channels make sound three-dimensional 
    f.setsampwidth(1) # Set the sampling bit width , That is, it takes up memory , For example, the bit width of lossless music is at least 2 byte , That is to say 0~65535
    f.setframerate(frameRate) # Tell the computer how many data a second 
    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)

Draw sound files ：

#!/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 Introduction to the algorithm ( Discrete Fourier transform )

effect ： Periodic signal decomposition （ Don Fourier thinks that all signals are periodic , So it can also be called signal decomposition ）

For example, a complex signal is generated by 1Hz and 2Hz Composed of ,DFT The algorithm is to decompose this complex signal into 1Hz and 2Hz

The formula ：
$$X(m)=\sum _{n=0}^{N-1}x(n)\cdot e^{-j2\Pi mn/N}$$

Formula analysis ：X(m) On behalf of the m Analysis results of three frequencies ,x（n） Medium n On behalf of the n Sampling signal data at time point , For the entire e Medium m representative m Two frequencies

$$\begin{gathered} repairchargeThe\; baseFoundationknowknowledge： \\ {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} have\; toTo： \\ X(m)=\sum _{N=0}^{N-1}x(n)\cdot (w_N^{-m}{)}^n \end{gathered}$$

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

FFT And DFT contrast ：

FFT Represents a method of data processing , It's an algorithm

DFT Represents a specific kind of problem , Used to decompose signals , Time domain data （ Data given in chronological order ） Convert to frequency domain data , For example, the matrix above x(n) It's time domain ,X(m) The frequency domain

actual combat 1： Sound restoration

Ideas ： Parsing sound data ->DFT Algorithm -> Analyze frequency domain data -> Restore sound

1、 Parsing sound data ：（ A string of digital files will be generated ）

#!/usr/bin/env python
# coding=utf-8

import wave
import numoty as np

f  = wave.open("xxx.wave", "rb");//xxx.wave For a complex sound file 
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 Algorithm ：

#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 representative sin The preceding symbol 
        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 ;
    }
    // Get a point value to represent 
    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];// preservation a
        __transform(temp, n, -1);// Follow fft Change to negative on the basis 
        return temp;
    }
    // Get a coefficient to represent 
    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);// stay fft Replace with positive on the basis 
        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、 Analyze frequency domain data

A diagram of the transformation of the result obtained by the above discrete Fourier transform

For the graph above , according to DFT The properties of the algorithm are symmetric , So only half of the graph is useful . Just look at the height of the column , The height of the column represents the size of the sound , The bottom is the frequency

4、 Restore sound

#!/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) # from 0 To start generating time It's worth , The difference between every two values is 1.0 / frameRate
fileName = "guess.wav"

#wave_data = volumn * np.sin(2.0 * pi * readlf * t);#readilf Is the true frequency domain ,volumn For sound size , These two are through DFT Got 
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) # Mono , Two channels make sound three-dimensional 
f.setsampwidth(1) # Set the sampling bit width , That is, it takes up memory , For example, the bit width of lossless music is at least 2 byte , That is to say 0~65535
f.setframerate(frameRate) # Tell the computer how many data a second 
f.writeframes(wave_data.tostring())
f.close();
print("write data to : " + fileName)

actual combat 2： Noise reduction

A complex sound results DFT Converted image :

reflection ： What frequency signal is noise ？ The sound is low and the frequency is high , Because the sound with too high frequency is inaudible to us , So it doesn't make any sense . So noise reduction is to filter out all high-frequency signals .

because DFT The result is symmetric , Therefore, the closer the above figure is to the middle point, the higher the high-frequency signal

expand ：FFT The algorithm can be used in signal processing and other scenarios , such as PS Cutout

The actual battle is not over yet