傅里叶变换

前言

\qquad 傅里叶变换，相信大家刚接触这个算法的时候会觉得比较难懂，我刚开始用的时候也是用不明白，后来翻阅教材、百度各种资料，终于明白了一点。最近闲来无事写给大家评评，望指正！
\qquad 这篇博客对初学者理解傅里叶变换的本质会有所帮助。文中一些数据并非严谨，若想深入研究请参考教材。

频率

\qquad 频率是个什么东西？如果细心观察你会发现频率无处不在。比如：一日三餐，你吃饭的频率是每天三次；每天太阳东升西落，每天太阳升起的频率是一；你每天起床、坐车、吃饭、工作，周而复始，这是你生活的频率。再比如：手机通信，卫星通信，雷达，生活中的电子产品，这一切都可以由频率去表示。
\qquad 我们可以用频率去表示很多东西，比如：近红外频率范围120THz～380THz；远红外频率范围300GHz～6000GHz；可见光中红色的频率范围385THz～482THz，还有近年来比较火的5G通信。我们可不可以反过来看呢？用频率去表示，120THz～380THz是近红外，300GHz～6000GHz是远红外，当然也不排除在这个频段内有其他的东西就像每个人可以重名一样。我们可以认为世间万物都有一个频率的标签（说的比较夸张^^），都可以由频率对它们表示。
\qquad 那么如何知道一个事物的频率是多少呢？比如：电路板上的晶振频率是多少？当然可以用示波器直接看周期就能计算出频率。如下图，看一下刻度就能算出频率。

但是下面这幅图中的频率有哪些呢？不知到了吧，傅里叶告诉你它知道！

傅里叶变换

傅里叶级数
\qquad 傅里叶变换源自于傅里叶级数展开，傅里叶级数告诉我们任何一个周期函数都可以展开成傅里叶级数，因此各种波形的周期信号都可分解为一系列不同频率的正弦波。这就告诉我们信号可以由不同频率的正弦波合成，当然也可以分解为不同频率的正弦波。
比如：方波如何合成？方波公式如下：
$$f(x)=\frac{4h}{pi}(sin(wt)+\frac{1}{3}sin(3wt)+\frac{1}{5}sin(5wt)+...+\frac{1}{2n-1}sin(2n-1)wt)$$
其中h代表幅值，我们用matlab合成看一下。取h=4，n=5；

上图可以看出是方波的样子，但还是比方波差一点，这是因为n只取到了5，如果n取无限大它就会无限接近方波，这里就不一一例举了。

傅里叶变换
由傅里叶级数推导傅里叶变换的过程这里就不叙述了，直接看傅里叶变换公式：
$$F(w)={\int }_{-\infty }^{+\infty }f(t)\ast e^{-jwt}dt$$
这里为了方便叙述就用了连续傅里叶变换公式。
傅里叶变换是将时域信号转到频域！时域、频域是个什么东西？

什么是时域？
\qquad 百度的解释：时域是描述物理信号对时间的关系。这里说物理信号有些片面。举个俗一点的例子，汽车行驶的距离与时间的关系，人的衰老程度与时间的关系，电器故障的概率与时间的关系，等等。其本质是以时间为参考去观察事物的特点。在坐标系中的表示是：以汽车行驶距离为例，横轴是时间，纵轴是距离。

什么是频域？
\qquad 与时域的解释类似。百度的解释：频域是描述信号在频率方面特性时用到的一种坐标系。在电子学，控制系统工程和统计学中，频域图显示了在一个频率范围内每个给定频带内的信号量。暴力一点的解释就是事物的特性与频率之间的关系，在坐标系中横坐标是频率，纵坐标是幅值。以人为例，人体红外线的频率大概是31.9THz（不是严谨的数值），如下图人体红外线视图。

那么在坐标系中表示就是：在横轴31.9THz的位置，纵轴幅值为非零。

什么是傅里叶变换呢？
\qquad 傅里叶变换就是时域信号乘上傅里叶变换因子转化到频域。其实我也不知道怎么解释，它就是这么简单的逻辑，就像：
\qquad\qquad 克拉克穿上红内裤就变成超人，脱下就变成了记者；
\qquad\qquad 彼得帕克穿上战斗服就是蜘蛛侠，脱下就是中学生；
\qquad\qquad 傅里叶变换也一样，时域信号“穿上”变换因子就变成了频域，“脱下”变换因子就是时域。
\qquad 它的过程只是一个数学计算，就像s=v*t，将时间与距离做了一次映射，傅里叶变换也是一样将频率w与幅值F做了一次映射。

傅里叶变换是怎么求频率的？
\qquad 举个简单的例子，我们有十把钥匙去开一个门，找出能打开门的钥匙，开门方式是每把钥匙分别去试开，如果打开了就认为钥匙与门能对应上，并且记录，再用下一把钥匙试开。

\qquad 傅里叶变换的过程与这类似，将时域信号$f(t)$看成门，频率$w_n$看成是钥匙，用不同的频率$w_n$与信号中频率配对（因为信号可看成由若干个不同频率的信号叠加成的），如果信号中有能够与频率$w_n$配对（频率相同）的频率则输出非零，如果没有能够配对的频率则输出零（或趋近于零）。当输出非零时就认为信号中有对应的频率$w_n$，这就是一个遍历过程。

傅里叶变换为什么能求出频率？如何求出频谱？
\qquad 上面叙述了傅里叶变换的过程，这里说明傅里叶变换是怎么求频率的。重写公式：
$$F(w)={\int }_{-\infty }^{+\infty }f(t)\ast e^{-jwt}dt$$
根据傅里叶变换求频率的过程，先假设w值然后进行“配对”。
1）当$w=w_0$时，由傅里叶公式得，求积分。其中f(t)为时域采集的信号，$T_p$为采样周期。
$$y(w_0)=F(w_0)={\int }_t^{t+T_p}f(t)\ast e^{-j{w}_{0}t}dt$$
2）当$w=w_1$时，同上
$$y(w_1)=F(w_1)={\int }_t^{t+T_p}f(t)\ast e^{-j{w}_{1}t}dt$$
3）当$w=w_n$时，同上
$$y(w_n)=F(w_n)={\int }_t^{t+T_p}f(t)\ast e^{-j{w}_{n}t}dt$$
通过上述过程，我们得到了一组数列$[(w_0,y(w_0))，(w_1,y(w_1))，(w_2,y(w_2))，\dots ，(w_n,y(w_n))]$
将$w_n$看成坐标横轴，$y(w_n)$看成坐标纵轴，就得到频谱图。它是一个对称谱。

为什么当$f(t)$中包含$w_0$频率时得到一个非零值呢？
现在需要了解几个公式：
欧拉公式
$$e^{jx}=cos(x)+jsin(x)$$
积化和差

$$\sin \alpha \cos \beta =\frac{1}{2}[\sin (\alpha +\beta )+\sin (\alpha -\beta )]$$
$$\cos \alpha \sin \beta =\frac{1}{2}[\sin (\alpha +\beta )-\sin (\alpha -\beta )]$$
$$\cos \alpha \cos \beta =\frac{1}{2}[\cos (\alpha +\beta )+\cos (\alpha -\beta )]$$
$$\sin \alpha \sin \beta =-\frac{1}{2}[\cos (\alpha +\beta )-\cos (\alpha -\beta )]$$

对傅里叶变换公式进行分解，假设信号$f(t)=Acos(w_{c}t)$，A为信号幅值，$T_p$为采样周期（为方便计算设$T_p$为pi的整数倍周期），则：

$$\begin{array}{rl}F(w_0)& ={\int }_t^{t+T_p}f(t)\ast e^{-j{w}_{0}t}dt\\ & =A{\int }_t^{t+T_p}\cos (w_{c}t)\ast e^{-j{w}_{0}t}dt带入信号\\ & =A{\int }_t^{t+T_p}\cos (w_{c}t)\ast (\cos (w_{0}t)-j\sin (w_{0}t))dt欧拉公式\\ & =A{\int }_t^{t+T_p}\cos (w_{c}t)\ast \cos (w_{0}t)dt-Aj{\int }_t^{t+T_p}\cos (w_{c}t)\ast \sin (w_{0}t)dt\\ & =\frac{1}{2}A{\int }_t^{t+T_p}(\cos (w_{c}t+w_{0}t)+\cos (w_{c}t-w_{0}t))dt\\ & -\frac{1}{2}Aj{\int }_t^{t+T_p}(\sin (w_{c}t+w_{0}t)-\sin (w_{c}t-w_{0}t))dt\\ 当w_c=w_0时& \\ & =\frac{1}{2}A{\int }_t^{t+T_p}(\cos (2w_{0}t)+\cos (0))dt-\frac{1}{2}Aj{\int }_t^{t+T_p}(\sin (2w_{0}t)-\sin (0))dt\\ & =0+\frac{1}{2}A\ast T_p-0-0正余弦信号整数倍周期内积分为0\\ & =\frac{1}{2}A\ast T_p\\ 结果取归一化为\frac{1}{2}A，即信号幅值的一半。& \end{array}$$

从推导结果可以看出当$w_c=w_0$时，即$f(t)$中拥有频率$w_0$时求得的结果不为零。
下面通过matlab证明：

clear all
close all
clc;

t=0.01:0.01:1;  %%采样点数100，采样频率为100Hz
Fc= 1;          %%信号频率为1hz  满足乃奎斯特定理
w0 = 2*pi*Fc;   %%角频率
A = 4;          %%幅值为4
ft=A*cos(w0.*t);%%信号模型
e0= exp(-1i*w0.*t);  %%傅里叶变换旋转因子
yt=sum(ft.*e0)/100;  %%积分运算，求归一化
NFFT = abs(fft(ft))/100;  %%fft运算,对比，取归一化
plot(NFFT);
xlim([-1 101]);
ylim([-0.1 2.2]);

可以看出结果实部是2，虚部是0 (10^(-17)可以看成是0)，求模的结果是2，这与理论计算吻合。

调用系统函数fft求得的结果是在第2个频点处信号幅值为2，这与上述结果吻合。

$$\begin{array}{rl}当w_c\ne w_0时& \\ 原式& =\frac{1}{2}A{\int }_t^{t+T_p}(\cos (w_{c}t+w_{0}t)+\cos (w_{c}t-w_{0}t))dt\\ & -\frac{1}{2}Aj{\int }_t^{t+T_p}(\sin (w_{c}t+w_{0}t)-\sin (w_{c}t-w_{0}t))dt\\ & =\frac{1}{2}A{\int }_t^{t+T_p}\cos ((w_c+w_0)t)dt+\frac{1}{2}A{\int }_t^{t+T_p}\cos ((w_c-w_0)t)dt\\ & -\frac{1}{2}Aj{\int }_t^{t+T_p}\sin ((w_c+w_0)t)dt+\frac{1}{2}Aj{\int }_t^{t+T_p}\sin ((w_c+w_0)t)dt\\ & =0\end{array}$$
下面通过matlab证明：

t1=0.01:0.01:1;         %%采样点数100，采样频率为100Hz
Fc= 1;                  %%信号频率为1hz  满足乃奎斯特定理
F1= 3;
wc = 2*pi*Fc;           %%信号角频率
w1 = 2*pi*F1;           %%旋转因子角频率
A = 4;                  %%幅值为4
fc=A*cos(wc.*t1);       %%信号模型
e1= exp(-1i*w1.*t1);    %%傅里叶变换旋转因子
yt1=sum(fc.*e1)/100;    %%积分运算，求归一化

结果是个非常小的值，可以认为是0。与上述推论吻合。
\qquad 上述过程中我取的采样时间t为pi的整数倍，有人会问为什么是整数倍？不是整数倍对不对呢？不是整数倍求得的结果也是对的，只是sin、cos求积分的时候不为0，然后就很难解释了，我为了避免给自己挖坑所以避开了这点^^。感兴趣的话自己用matlab算一下吧。

\qquad 在上述分析中为了方便用的信号是余弦信号，实际工程中用的信号为复指数信号$e^{jwt}$，其实复指数信号跟正余弦信号是一样的，根据欧拉公式复指数信号可以转化为正余弦信号。

\qquad 上述分析就是此篇博客的主要内容，下面再多少几句废话^^

离散傅里叶变换

\qquad 上面的叙述应用的都是连续傅里叶变换，在实际应用中用的都是离散傅里叶变换，那么什么是离散傅里叶变换呢？
\qquad 离散傅里叶变换可以认为是连续傅里叶变换的频域采样。下面是离散傅里叶变换的基本公式：
$$F(w)=\sum _{n=-\infty }^{+\infty }x(n)e^{-jwn}$$
对公式进行变形，取
$$w=2\ast pi\ast f$$
取采样点数为N，$f^{\prime }=\frac{1}{N}$
$$w_T=2\ast pi\ast f^{\prime }=2\ast pi\ast \frac{1}{N}$$
上述的$w_T$是离散傅里叶变换的基频，即每次取的频率都是$w_T$的整数倍。
$$w=w_T\ast k=2\ast pi\ast \frac{1}{N}\ast k$$
$$F(w_T\ast k)=\sum _{n=-\infty }^{+\infty }x(n)e^{-j\ast 2\ast pi\ast \frac{1}{N}\ast k\ast n}$$
也可以写成
$$F(k)=\sum _{n=-\infty }^{+\infty }x(n)e^{-j\ast 2\ast pi\ast \frac{1}{N}\ast k\ast n}$$
这样求出的结果可以说成是第k个频点的值是多少。其频谱分辨率是$w_T$。

FFT是个啥东西？

\qquad 从前面的叙述可以发现傅里叶变换是个遍历过程，计算量非常大。在傅里叶变换刚发现的时候由于计算量太大而无法应用，后来通过研究它的性质发现的一种快速算法使计算量大大降低，这就是快速傅里叶变换FFT。FFT算法这里就不细说了。

结束

\qquad 以上就是我认为的傅里叶变换，其实它没有那么复杂，它就是一个公式，带入参数计算得到结果，就这么简单。
推荐一本书吧，<<数字信号处理—原理、算法与应用>>（第四版）作者 John G.Proakis Dimitris G.Manolakis。里面有详细的理论推导，想要研究的可以看看。