【今日话题】在嘈声中精确分离出个体——“正所谓:长乐钟声花外尽,龙池柳色雨中深

2023-12-16 16:43:54

傅里叶变换

系数表示法与点值表示法

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)0n?f(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] ?y0?y1?y2?y3?...yn?? ?= ?11111?x0?x1?x2?x3?xn??x02?x12?x22?x32?xn2??x03?x13?x23?x33?...xn3?????????,x0m??,x1m??,x2m??,x3m??,xnm?? ?× ?a0?a1?a2?a3?...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/2 ?i?1/2 ?+1/2 ?i?ii?1/2 ??1/2 ?i1/2 ?+1/2 ?i?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/2 ?i?1/2 ?+1/2 ?i?ii?1/2 ??1/2 ?i1/2 ?+1/2 ?i?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/2 ?i
和与其对称的点
? 1 / 2 ? 1 / 2 i -\sqrt{1/2}-\sqrt{1/2}i\\ ?1/2 ??1/2 ?i

在这里插入图片描述

同理对-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/2 ?i?1/2 ?+1/2 ?i

总结:要求某个点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(16?)y=sin(16?)

那么就可以得到公式:(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(nk?)+sin(nk?)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(nk?)?sin(nk?)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?)+wnk?P1?(wn/2k?)f(wnk+n/2?)=P0?(wn/2k?)?wnk?P1?(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?)+wnk?P1?(wn/2k?)f(wnk+n/2?)=P0?(wn/2k?)?wnk?P1?(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] ?P0?P1?P2?P3?? ?= ?1111?w40?w41?w42?w43??(w40?)2(w41?)2(w42?)2(w43?)2?(w40?)3(w41?)3(w42?)3(w43?)3? ? ?C0?C1?C2?C3?? ?

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/4?w40?/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? ? ?P0?P1?P2?P3?? ?= ?C0?C1?C2?C3?? ?

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? ?1111?w40?w4?1?w4?2?w4?3??(w40?)2(w4?1?)2(w4?2?)2(w4?3?)2?(w40?)3(w4?1?)3(w4?2?)3(w4?3?)3? ? ?P0?P1?P2?P3?? ?= ?C0?C1?C2?C3?? ?

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

若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=0N?1?x(n)?e?jmn/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 补充基础知识:ejk/n=wnk?=cos(nk?)+sin(nk?)ie?jk/n=wn?k?=cos(nk?)?sin(nk?)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=0N?1?x(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] ?1111?w40?w4?1?w4?2?w4?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抠图

实战未完待续
在这里插入图片描述

文章来源:https://blog.csdn.net/qq_58360406/article/details/135033220
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。