从傅里叶级数到基于互相关运算的声音测距【学习笔记】

傅里叶变换

傅里叶级数

任何周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示
设一函数周期为$2l$
$$f(x)=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_{n}cos(\frac{n2\pi x}{2l})+\sum _{n=1}^{\infty }{b}_{n}sin(\frac{n2\pi x}{2l}),n\in N^{\ast }$$
可以找到一系列特定周期的正余弦函数，使它们互相正交，即乘积在最大周期内的积分为零
${\int }_{-l}^{l}sin(\frac{n\pi x}{l})dx=0;$
${\int }_{-l}^{l}sin(\frac{m\pi x}{l})cos(\frac{n\pi x}{l})dx=0,m,n\in N^{\ast };$
${\int }_{-l}^{l}sin(\frac{m\pi x}{l})sin(\frac{n\pi x}{l})dx=0,m,n\in N^{\ast };$
……
据此可以以{$1,sin(\frac{n\pi x}{l}),cos(\frac{n\pi x}{l})$}为基表示出一个周期函数f(x)
求f(x)与每一个基的内积，即相乘求积分，即可得到每个基所对应的系数
$$a_0=\frac{1}{2l}{\int }_{-l}^{l}f(x)dx,$$
$$a_n=\frac{1}{l}{\int }_{-l}^{l}f(x)sin(\frac{n\pi x}{l})dx,$$
$$b_n=\frac{1}{l}{\int }_{-l}^{l}f(x)cos(\frac{n\pi x}{l})dx,$$
$\frac{n\pi }{l}=2\pi f$,若将频率$f$作为自变量，将$a_n,b_n$作为因变量，便可以得到一个频域的图像，有助于我们分析一段信号里的频率特征

下面引入欧拉公式，可以将傅里叶级数写成复数形式：
$$f(x)=\sum _{n=-\infty }^{+\infty }{c}_n{e}^{i\frac{n\pi x}{l}}$$
$$c_n=\frac{1}{2l}{\int }_{-l}^{l}f(x)e^{-i\frac{n\pi x}{l}}dx$$

对于非周期函数，直接以周期无穷大处理，而对于定义域连续且有限的一段信号，可以以定义域长度为周期，进行傅里叶变换

---

离散傅里叶变换(DFT)

现实生活中，对信号的数字化采样是离散的，这时就需要离散形式的傅里叶变换（DFT）
主要的改变就是把积分形式改为求和
设有一信号采样频率为$f_s$,采样点数为N，那么对它进行傅里叶变换后可以得到N个频率幅值
$$c_n=\frac{1}{N}\sum _{k=0}^{N-1}f(k)e^{i\frac{2\pi n}{N}k}$$

• 通常也把$e^{i\frac{2\pi }{n}k}$记作${\omega }_n^k$ , 之后推导FFT的时候会用到
• （余以为，指数上的负号应该是用欧拉公式推导复数形式时产生的，这里直接用实部记录余弦分量，以虚部记录记录正弦分量，便省去负号）

对 n ∈ { 0 , 1 , . . . , N − 1 } n\in\{0,1,...,N-1\} n∈{0,1,...,N−1}依次计算，对得到的$c_n$取模， 便可以得到一个频率从0~$f_s$的频谱图。自变量频率是等间距分布的，即$c_n$对应的实际频率为 f s ∗ n N , n ∈ { 0 , 1 , . . . , N − 1 } f_s*{n\over N},n\in\{0,1,...,N-1\} fs​∗Nn​,n∈{0,1,...,N−1}
同时采样定理告诉我们若要完整地保留原始信号中的信息，采样频率应至少大于原始信号中最高频率的两倍。可以只采用前一半结果。（具体原因不懂我瞎说的）
综上，采样频率越高，可以准确测量的频率也就越高（采样频率的一半）；参与DFT运算的采样点越多，频谱的分辨率（$\frac{N}{f_s}$）就越高，同时也必然导致更长的采样时间和运算时间。

至此已经可以写出离散傅里叶变换(DFT)的代码：

#define N 1024
#define PI 3.1415926535
double x[N];//输入信号 
double Re[N],Im[N];//实部虚部分别为cos、sin
double c[N];//输出频谱
int n,k;
for(n=0;n<N;n++)
{
   for(k=0;k<N;k++)
   {
   	Re[n] += x[k]*cos(2*PI*k*n/N);
   	Im[n] += x[k]*sin(2*PI*k*n/N);
   }
   c[n] = sqrt(Re[n]*Re[n]+Im[n]*Im[n]);//取模
}

这种写法十分简洁易懂，但是很费时。在I.MX RT1064上测试计算1024个采样点需要约1000ms

---

快速傅里叶变换(FFT)

快速傅里叶变换是一种可以加速离散傅里叶(DFT)运算速度的算法，大大减少了运算量。同时得到傅里叶变换的准确结果。

单位根及其性质

单位根
在复平面的单位圆上取n个点，这些点将圆周n等分。那么它们的指数表示为 e i 2 π n k ， k ∈ { 1 , 2 , . . . , n − 1 } e^{i{2\pi \over n}k}，k\in \{1,2,...,n-1\} ein2π​k，k∈{1,2,...,n−1}，现将它们记为 ω n k ， k ∈ { 1 , 2 , . . . , n − 1 } \omega_n^k，k\in \{1,2,...,n-1\} ωnk​，k∈{1,2,...,n−1}，直观来看下标就是圆周被等分的份数，上标是从实轴开始逆时针数过的圆弧个数。

显然，这些复数的模长为1，且它们的n次方都为实数1，故称之为单位根。

$${\omega }_n^k=e^{i\frac{2\pi }{n}k}$$

单位根的性质

$${\omega }_n^k={\omega }_{2n}^{2k}$$

将圆弧再等分依次，变成2n等分，取2k份圆弧，两者表示同一个复数。也可以直接通过指数形式理解。

$${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$$

$+\frac{n}{2}$表示在单位圆上逆时针旋转了半圈，即与原点对称，实部虚部取相反数。

$$({\omega }_n^k{)}^m={\omega }_n^{mk}$$

这可以由复数乘法性质得到，模长相乘(=1),幅角相加。

---

单位根如何加速DFT运算

再贴一下离散傅里叶的公式：
$$c_n=\frac{1}{N}\sum _{k=0}^{N-1}f(k)e^{-i\frac{2\pi n}{N}k}$$
将$e^{i\frac{2\pi }{n}k}$记为${\omega }_n^k$形式：
$$c_n=\frac{1}{N}\sum _{k=0}^{N-1}f(k){\omega }_N^{nk}$$
或$$c_n=\frac{1}{N}\sum _{k=0}^{N-1}f(k)({\omega }_N^n{)}^k$$
这是我们需要计算的任务。
设$A(x)={\sum }_{k=0}^{N-1}f(k)x^k，$并将$f(k)用a_k$表示
则$$A(x)=\sum _{k=0}^{n-1}{a}_k{x}^k=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$$
这里开始采用分治的思想，按下标奇偶性将数列分成两半（此处规定n为偶数）：
$$\begin{gathered} A(x)=(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+...+a_{n-1}{x}^{n-1}) \\ =(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+...+a_{n-1}{x}^{n-2}) \end{gathered}$$
令$$\begin{gathered} A_1(x)=(a_0+a_{2}x+...+a_{n-2}{x}^{\frac{n-2}{2}}), \\ {A}_2(x)=(a_1+a_{3}x+...+a_{n-1}{x}^{\frac{n-2}{2}}). \end{gathered}$$
则$$A(x)=A_1(x^2)+x{A}_2(x^2)$$
接下来开始带入单位根${\omega }_n^k$(这里令$k<\frac{n}{2}$,分治减少运算量的关键步骤):
$$\begin{gathered} A({\omega }_n^k)=A_1(({\omega }_n^k{)}^2)+{\omega }_n^k{A}_2(({\omega }_n^k{)}^2) \\ =A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k}) \\ =A_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$
对于$k>\frac{n}{2}$，不妨直接带入${\omega }_n^{k+\frac{n}{2}}$(多用单位根的几何意义理解):
$$\begin{gathered} A({\omega }_n^{k+\frac{n}{2}})=A_1(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}{A}_2(({\omega }_n^{k+\frac{n}{2}}{)}^2) \\ =A_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2k+n}) \\ =A_1({\omega }_n^{2k}\ast {\omega }_n^n)-{\omega }_n^k{A}_2({\omega }_n^{2k}\ast {\omega }_n^n) \\ =A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k}) \\ =A_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$
观察以上结果，发现前一半($k<\frac{n}{2}$)和后一半的结果十分相似，仅仅差一个符号
根据这个结论，我们就可以只带入前一半单位根的同时得到整个序列的DFT结果

再看两个推导的最后一行，$n写作了\frac{n}{2}$,有了这一步，就可以发现$A_1({\omega }_{\frac{n}{2}}^k)、A_2({\omega }_{\frac{n}{2}}^k)和A({\omega }_{\frac{n}{2}}^k)$其实是等价的，只不过取的单位根少了一半，而这两个一半长度的序列又可以通过上面的分治的方法减半运算……
（典型的递归思想，这时候就需要满足$n=2^i,i\in N^{\ast }$了）

总得来说，FFT是充分利用了复数表示中复数的特性（见单位根的性质），推导出了分治的计算方法，并且可以递归使用，从而大大减少了计算复杂度

---

递归实现
现在就可以用递归的方式完成FFT的代码：
C语言的复数库对于不同编译器使用不太相同，所以先写一个简单的复数计算需要用到的代码

typedef struct complex_t
{
    double x,y;
} complex;

complex plus(complex a,complex b)//+
{
    a.x += b.x;
    a.y += b.y;
    return a;
}

complex minus(complex a,complex b)//-
{
    a.x -= b.x;
    a.y -= b.y;
    return a;
}

complex multi(complex a,complex b)//*
{
    complex m;
    m.x = a.x*b.x - a.y*b.y;
    m.y = a.x*b.y + a.y*b.x;
    return m;
}

接下来是递归程序

#define N 1024  //序列总长度，需满足2的整数次幂
#define Pi 3.14159265359

void FFT_recursion(int len,complex *x)//当前序列长度、当前求解序列
{
    int i;
    complex x1[N],x2[N];	//放按下标奇偶分组的两个序列
    complex Wn = {cos(2.0*Pi/len),sin(2.0*Pi/len)};	//k=1单位根，计算A2前系数用
    complex w = {1,0},m;	//w用来保存单位根的幂
    if(len == 1)return;	//当长度为1时就可以直接返回，什么也不做
    for(i=0;i<len;i+=2)		//给半长序列填入相应的值
    {
         x1[i>>1] = x[i];
         x2[i>>1] = x[i+1];
    }
    FFT_recursion(len>>1,x1);	//求出两个序列分别带入len/2个单位根的值
    FFT_recursion(len>>1,x2);	//结果就存储在输入的数组中
    
    for(i=0;i<(len>>1);i++,w = multi(w,Wn))	//这里i相当于公式中的k
    {
         m = multi(w,x2[i]);	//保存w与A2的乘积
         x[i] = plus(x1[i],m);	//同时给前后两部分赋值
         x[i+(len>>1)] = minus(x1[i],m);
    }
}

如果使用递归写法，就需在调用函数后再对类型为complex的数组进行处理
可以看到这里求的是$A({\omega }_n^k)$，要求严谨的幅值还需要再除以N，并且取模方便使用（其实工程中直接除以一个合适的值，方便观察就可以了）

用递归实现的FFT已经比之前的DFT高效很多了，但需要特殊功能的时候还得用另外的函数封装，大数组又占空间，如果能改成用for循环实现岂不美哉？

---

迭代实现
用递归函数写代码的时候，我们只需要思考其中一步需要完成的操作，以及回归的条件。
那么最终整个递归过程都做了些什么呢，下面我们以N=8为例，把每一步分治的过程写出来。

第一行是原始输入的一组数据，我们要做的是用它们分别乘上${\omega }_8^k,k\in [0,7]$然后求和
第二行采用分治方法，按奇偶性分成两列，对它们求值后得到两个长度4的数组，设为A1[ ],A2[ ]，那么用A1[0]±$({\omega }_8^1{)}^0$就可以得到最终需要计算的A[0]和A[4]两个值，其余三对结果求法类似。

下面把第一行和最后一行的下标及其二进制写出来：

原序列	0	1	2	3	4	5	6	7
二进制	000	001	010	011	100	101	110	111
重排后	0	4	2	6	1	5	3	7
二进制	000	100	010	110	001	101	011	111

可以看出排列后递归最底层的顺序是有规律的——下标的二进制表示颠倒

再放一张N=16的表格：

原序列	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
二进制	0000	0001	0010	0011	0100	0101	0110	0111	1000	1001	1010	1011	1100	1101	1110	1111
重排后	0	8	4	12	2	10	6	14	1	9	5	13	3	11	7	15
二进制	0000	1000	0100	1100	0010	1010	0110	1110	0001	1001	0101	1101	0011	1011	0111	1111

下面研究一下如何用程序完成这个排列
（不想研究的可以抄代码，反正就一行 ）
如果把下标的二进制当成数组，一位一位地交换，随着n的增大会越来越麻烦

• c语言中右移可以表示整除2，比如5>>1=2, 6>>1=3
• 一个数（除了二进制最低位）可以看作它整除2的数左移一位，如13(1101)，而6(0110)<<1=12(1100)，前三位和13相同
• 一个二进制数颠倒后（除了最高位），可以看作它整除2的数也颠倒后右移一位，如13(1101) -颠倒-> 11(1011), 6(0110) -颠倒-> 6(0110) >> 1 = 3(0011)，后三位相同
• 剩下的一位单独取出（和1作逻辑与运算）进行位移后和其它位作或运算即可
• 我们只要从0(0本身就对称，移位后也无变化)开始依次操作，每次去取整除2（右移1）的下标所对应的排列后标号，就能把所有排列计算出来

下面是代码：

#define N 1024	//总数据个数
#define L 10	//满足 N = 2^L	//同时也是下标二进制表示的位数

int r[N] = {0};	//可以认为下标是原排列，数据是二进制颠倒后的排列
int i;	//循环变量

for(i=0;i<N;i++)
{
    r[i] = (r[i>>1]>>1)|(i&1)<<(L-1); //计算重排后的数组下标
    //逻辑或符号左侧是取下标整除2的结果右移一位
    //右侧是单独处理剩下来的一位，将下标本身第一位取出，左移L-1位
}
//这部分代码可以直接放进迭代FFT函数中

现在我们已经成功得到了一个排列后的下标，只要按照这个数组给出的顺序循环填入数据就可以开始进行计算了

直接给出完整代码吧：

#define N 1024 //总数据个数
#define L 10 //满足 N = 2^L
#define Pi 3.14159265359

void FFT(complex *x)
{
    int r[N] = {0},i,l,j,k;
    complex Wn,w;
    complex temp,temp1;	//用于交换数组项，以及计算时保存数组原值
    for(i=0;i<N;i++)
        r[i] = (r[i>>1]>>1)|(i&1)<<(L-1); //计算重排后的数组下标

    for(i=0;i<N;i++)	//按r[i]所给顺序重排数组
        if(i<r[i])	//防止重复交换
        {
            temp = x[i];
            x[i] = x[r[i]];
            x[r[i]] = temp;
        }
    
    for(l=1;l<N;l<<=1)	//当前需要计算的数组长度的【一半】，l取1,2,4,8,...,N/2
    {
        Wn.x = cos(Pi/l);Wn.y = sin(Pi/l);	//准备单位根
        for(j=0;j<N;j+=l<<1)	//依次处理每个序列，j每次循环递增2l
        {
            w.x = 1;w.y = 0;	//记录单位根的幂
            for(k=0;k<l;k++,w = multi(w, Wn))	//用上一步算好的值带入公式计算当前序列
            {
                //由于迭代法只用了一个复数数组，结果也直接记录在本身，
                //所以要先用临时变量记录原值。
                temp = x[j+k];temp1 = multi(w,x[j+k+l]);
                x[j+k] = plus(temp, temp1);
                x[j+k+l] = minus(temp, temp1);
            }
        }
    }
}

如需作频谱分析，测得的信号填入数组实部，输出数组进行取模（平方和开根号）即可
实测以上代码在I.MX RT1064上仅需1.5ms！

至此FFT的研究就基本完成了！

---

多项式乘法

逛了几天，看到最多的对FFT作用的描述就是：加速多项式乘法。
苦恼了好长时间，始终想不通傅里叶级数里哪里用到多项式了。最后才明白，其实是用傅里叶变换的快速算法，去加速多项式的乘法计算。这里已经不谈FFT频域转换的功能了。

多项式的表示

讲多项式乘法还是离不开两种表示方法

系数表示
下面是一个普通的多项式，它有n项，次数是n-1
$$f(x)=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$$
把它的系数拿出来，列成一个向量，就得到了多项式的系数表示：
$$(a_0,a_1,a_2,...,a_{n-1})$$
这是符合我们对多项式的一般认识的表示方法

点值表示
我们选一些x的值带入表达式中，就可以求得多项式的值。如果带入了n个x，得到n个值，理论上就可以通过这一系列方程求解出原式多项式的每一个系数。
设 $yk=f(x_k)$，然后带入 x k , k ∈ { 0 , 1 , . . . , n − 1 } x_k,k\in \{0,1,...,n-1\} xk​,k∈{0,1,...,n−1}
得$$(y_0,y_1,y_2,...,y_{n-1})$$
这就是点值表示法
点值表示的多项式作乘法运算时只需要n次乘法。

---

FFT如何加速多项式乘法

FFT作多项式乘法的过程主要包括求值(DFT)，插值(IDFT) 两部分。
求值指带入x的值，将系数表示转为点值表示；插值指将点值表示转回系数表示。

所以，求两个n项n-1次多项式相乘的结果，首先要对它们分别带入n个x值，求出点值表达（求值），然后循环n次将点值表示相乘，最后把相乘的结果还原为系数表达（插值）。

听起来是不是十分复杂，甚至有点多此一举？而且插值过程可是要求解n元一次方程组啊！
不过要知道，FFT对DFT的加速可不止一倍两倍，而且如果我们求值的过程中选取上述的单位复根作为自变量带入，插值过程就完全不需要解方程！而且与求值过程仅仅一个符号只差（我看呆了 ）！

---

求值

求值过程与之前讲的FFT类似，同样是带入单位复根${\omega }_n^k$，代码完全相同。
$$y[k]=\sum _{i=0}^{n-1}x[i]({\omega }_n^k{)}^i$$
x[ ]为原多项式的系数，y[ ]记录点值表示。而代码中结果直接输出到x[ ]本身。
这里直接调用该函数：

void FFT(complex *x);	//输入复数组，输出求值结果

---

插值

插值过程与点值过程的区别就在于：

1. 带入自变量时取单位复根的共轭复数 ${\omega }_n^{-k}$
2. 结果要除以n

插值过程所作的事情见以下表达式：
$$c_k=\frac{1}{n}\sum _{i=0}^{n-1}y[i]({\omega }_n^{-k}{)}^i$$
$c_k$为所求得的系数，$y[i]$为已经做过乘法的点值表达
可以看到这里似乎又把已求得的点值表达当成了系数带入求值了，为什么这样做就能解出系数呢？（又要敲公式了，不想看的可以跳过 ）

以下是证明：
设所求多项式的系数表达为$(a_0,a_1,a_2,...,a_{n-1})$，现有的是它经过FFT求值后的点值表达$(y_0,y_1,y_2,...,y_{n-1})$，（其实是通过前面求值的点值多项式相乘得到的）

现设一向量$(c_0,c_1,c_2,...,c_{n-1})$满足$c_k=\frac{1}{n}{\sum }_{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$，我们要证明$c_k=a_k$.

首先带入$y_i={\sum }_{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j$，
$$\begin{gathered} c_k=\frac{1}{n}\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i \\ =\frac{1}{n}\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i \\ =\frac{1}{n}\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i)({\omega }_n^{-k}{)}^i \\ =\frac{1}{n}\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i) \\ =\frac{1}{n}\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i) \end{gathered}$$
调换求和次序
$$c_k=\frac{1}{n}\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)$$
当$j=k$时，${\omega }_n^{j-k}=1$，$({\sum }_{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)=n$
当$j!=k$时，记${\omega }_n^{j-k}$为${\omega }_n^m$
$$\begin{gathered} \sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=1+({\omega }_n^m)+({\omega }_n^m{)}^2+...+({\omega }_n^m{)}^{n-1} \\ ({\omega }_n^m)\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=({\omega }_n^m)+({\omega }_n^m{)}^2+({\omega }_n^m{)}^3+...+({\omega }_n^m{)}^n \end{gathered}$$
做差，得
$$({\omega }_n^m-1)\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=({\omega }_n^m{)}^n-1$$
$$\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=\frac{({\omega }_n^m{)}^n-1}{{\omega }_n^m-1}$$
注意$({\omega }_n^m{)}^n=1$
所以${\sum }_{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=0$

再看上面的式子：
$$\begin{gathered} c_k=\frac{1}{n}\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i) \\ =\frac{1}{n}{a}_k(\sum _{i=0}^{n-1}({\omega }_n^0{)}^i) \\ =\frac{1}{n}{a}_k\ast n \\ =a_k \end{gathered}$$

---

证明完毕，可以开始写代码了！

#define N 8	//大于等于项数和且满足2的整数次幂
#define L 3

double* FFT_Convolution(double *a, double *b)
{
    int r[N] = {0},i,j,l,k;
    complex x[N],y[N];
    complex w,Wn,x1,x2,y1,y2;
    for(i=0;i<N;i++)
    {
        r[i] = (r[i>>1]>>1)|(i&1)<<(L-1); //计算重排后的数组下标
        x[i].x = a[r[i]];x[i].y = 0.0;	//按r[i]重排x,y
        y[i].x = b[r[i]];y[i].y = 0.0;
    }
    //数组次序排列完成
    for(l=1;l<N;l<<=1)  //求值
    {
        Wn.x = cos(Pi/l);Wn.y = sin(Pi/l);
        for(j=0;j<N;j+=l<<1)
        {
            w.x = 1.0;w.y = 0.0;
            for(k=0;k<l;k++,w = multi(w, Wn))
            {
                x1 = x[j+k];
                y1 = y[j+k];
                x2 = multi(w,x[j+k+l]);
                y2 = multi(w,y[j+k+l]);
                x[j+k] = plus(x1, x2);     //同时操作两个序列
                x[j+k+l] = minus(x1, x2);
                y[j+k] = plus(y1, y2);
                y[j+k+l] = minus(y1, y2);
            }
        }
    }
    
    for(i=0;i<N;i++)    //求积
        y[i] = multi(x[i],y[i]);
    
    for(i=0;i<N;i++)	//再次重拍
    	x[i] = y[r[i]];
    
    for(l=1;l<N;l<<=1)  //插值
    {
        Wn.x = cos(Pi/l);Wn.y = -sin(Pi/l);	//这里多了一个负号
        for(j=0;j<N;j+=l<<1)
        {
            w.x = 1.0;w.y = 0.0;
            for(k=0;k<l;k++,w = multi(w, Wn))
            {
                x1 = x[j+k];
                x2 = multi(w,x[j+k+l]);
                x[j+k] = plus(x1, x2);
                x[j+k+l] = minus(x1, x2);
            }
        }
    }
    
    for(i=0;i<N;i++)	//将结果填入数组b
        b[i] = x[i].x/N;

    return b;
}

• 需要注意的是N必须要大于等于所输入的两个多项式项数和

---

卷积与互相关运算

现在我们已经能求出多项式相乘的结果了，那么它有什么用呢？显然不仅仅是求多项式。

卷积

当我们平时手算多项式乘法时，如果次数比较高，求出来的项会特别多，所以我一般会按次数顺序计算，即先从0次项，再去找1次项的系数，以此类推。
如果把每一项的系数按顺序写出来，就会发现只要把两个系数相向放置，然后一项一项错开，把对应的系数相乘再求和就可以得到该次数的项的系数。
$$\begin{gathered} a_3,a_2,a_1,a_0.............. \\ ...............b_0,b_1,b_2,b_3 \end{gathered}$$$$a_0\ast b_0=c_0（常数项系数）$$
$$\begin{gathered} a_3,a_2,a_1,a_0.......... \\ .........b_0,b_1,b_2,b_3 \end{gathered}$$$$a_1\ast b_0+a_0\ast b_1=c_1（一次项系数）$$
$$\begin{gathered} a_3,a_2,a_1,a_0.... \\ .....b_0,b_1,b_2,b_3 \end{gathered}$$$$a_2\ast b_0+a_1\ast b_1+a_0\ast b_2=c_2（二次项系数）$$
$$......$$

即$$c_k=\sum _{i=0}^k{a}_i\ast b_{k-i}$$
这个操作就是卷积
卷积在工程信号处理中有着广泛的应用，遇到卷积时，我们就可以用上述FFT算法来计算。

double* FFT_Convolution(double *a, double *b);

---

互相关运算

与卷积的区别在于第一个序列不需要逆序
将序列 { a n } \{a_n\} {an​}固定，移动 { b n } \{b_n\} {bn​}，从-(n-1)到(n-1)

最直观的物理意义就是经过位移，两个序列中相似的部分重叠的时候，得到的运算结果最大。由此可以作信号的识别。

在MATLAB中的函数为xcorr(x,y);
在上述卷积代码中只需在数组b赋值时反向填入即可

下面是一个例子，可以明显看出第二个序列向左移动600个单位时和第一个序列最相似（这里是重叠）

---

声音测距

下图是一段Chirp信号

A chirp is a signal in which the frequency increases (up-chirp) or decreases (down-chirp) with time.
Chirp信号是一段频率随时间增加或降低的信号。

• 为什么用Chirp信号测距呢？
  声音测距的基本原理就是测量声音在空气中传播的延迟，通过发出一段Chirp信号，再去对它进行采集，通过与原信号进行互相关运算就可以知道它延迟的时间，从而算出距离。
  那为什么不直接判断相位差呢？可以简单的计算一下，常温中声速约为334m/s，假设声音信号传播10米后被接收，这期间过去了大约30ms，而一个440hz的声音信号已经走过了13个周期，这中间跳过的周期是我们无法测量的。
  而使用Chirp信号，我们可以按需设置它一个周期的时长，就可以避免这个问题，而且调频信号一个周期内的特征点比一个简单的正弦波要多得多，测量也就更加精确。

为了得到精确的时间延迟，也就是希望信号相关结果出现的峰值约尖锐越好，作为测距的声音信号需要频谱较宽，比如时间很短的脉冲信号、具有变频的Chirp信号、白噪声信号等等。利用声音定位的动物们常常使用的是Chirp信号。

采样完整周期的声音信号后，和原式数据作互相关运算，得到峰值处的横坐标，乘上采样周期即可得到声音传播延迟，再乘上声速即可算出距离。

---

终于终于结束了！（累die ）
这一周在各大论坛逛了无数篇有关傅里叶的博客，但奈何数学已经还了好多给老师，废了很大的劲才勉强理解。得出的结论是，每个人自身薄弱的地方不同，学习时着重的点也不同，写的东西大多是针对自身认为的难点。于是决定自己针对我理解中的难点，结合自己的理解写一篇笔记，以便日后查阅。（如果能帮到大家理解就更好啦）
其中一些内容包含自身理解，可能有错误或者不严谨的地方，查阅时要注意。

下面贴上主要参考的博文：

• 傅里叶级数 三角形式 到 复数形式
• 离散傅里叶变换-DFT（FFT基础）
• FFT频谱分析原理
• FFT算法讲解——麻麻我终于会FFT了！
• 快速傅里叶变换(FFT)详解
• 十分简明易懂的FFT（快速傅里叶变换）

——2020/3/25