离散傅⾥叶变换（DFT）
  对于第⼀幅图来说，它侧重展⽰傅⾥叶变换的本质之⼀：叠加性，每个圆代表⼀个谐波分量。第⼆幅图直观的表⽰了⼀个周期信号在时域与频域的分解。傅里叶变换公式证明
周期信号的三⾓函数表⽰
  周期信号是每隔⼀定时间间隔，按相同规律⽆始⽆终重复变化的信号。任何周期函数在满⾜狄利克雷条件下（连续或只有有限个间断点，且都是第⼀类间断点；只有有限个极值点），都可以展开成⼀组正交函数的⽆穷级数之和。使⽤三⾓函数集的周期函数展开就是傅⾥叶级数。对于周期为T 的信号f(t)，可以⽤三⾓函数集的线性组合来表⽰，即
f(t)=a_0+\sum_{n=1}^{\infty }(a_n\cos n\omega t+b_n\sin n \omega t)
  式中\omega=\frac{2\pi}{T}是周期信号的⾓频率，也成基波频率，n\omega称为n次谐波频率；a_0为信号的直流分量，a_n和b_n分别是余弦分量和正弦分量幅度。根据级数理论，傅⾥叶系数a_0、a_n、b_n的计算公式为：
\left\{\begin{matrix}a_0=\frac{1}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)dt \\ a_n=\frac{2}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)\cos{n\omega
t}dt,n=1,2,3,... \\ b_n=\frac{2}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)\sin{n\omega t}dt,n=1,2,3,... \end{matrix}\right.
  若将式⼦中同频率的正弦项和余弦项合并，得到另⼀种形式的周期信号的傅⾥叶级数，即
f(t)=A_0+\sum_{n=1}^{\infty}A_n\cos(n\omega t+\varphi_n)
  其中，A_0为信号的直流分量；A_1\cos(\omega t+\varphi_1)为信号的基频分量，简称基波；A_n\cos(n\omega t+\varphi_n)为信号的n次谐波，n ⽐较⼤的谐波，称为⾼次谐波。上式说明任何周期信号只要满⾜狄利克雷条件，都可以分解成直流分量和⼀系列谐波分量之和，这些谐波分量的频率是周期信号基频的整数倍。
  ⽐较两种三⾓函数形式的傅⾥叶级数，可以看出它们的系数有如下关系：
\left\{\begin{matrix}A_0=a_0 \\A_n=\sqrt{a_n^2+b_n^2},\quad\varphi_n=-\arctan\frac{b_n}{a_n} \\a_n=A_n\cos \varphi_n,\quad
b_n=A_n\sin\varphi_n \end{matrix}\right.
周期信号的复指数表⽰
  利⽤欧拉公式
\cos n\omega t=\frac{e^{jn\omega t}+e^{-jn\omega t}}{2}\quad\quad\sin n\omega t=\frac{e^{jn\omega t}-e^{-jn\omega t}}{2j}
  可以得到周期信号的复指数形式的傅⾥叶级数展开，即
  f(t)=\sum_{n=-\infty}^{\infty}F_ne^{jn\omega t}
  其中展开系数为
  F_n=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-jn\omega t}dt
  虽然n的取值范围为(-\infty,\infty)，但n取负值并不表⽰实际上存在负频率。周期信号可以⽤三⾓函数形式的傅⾥叶级数表⽰，也可以⽤复指数形式的傅⾥叶级数表⽰。前者物理意义明确，后者在数学处理上简便，并且容易与傅⾥叶变换统⼀起来。
周期信号的频谱图
  为了⽅便和直观的表⽰⼀个周期信号中所含有的频率分量，常⽤周期信号各次谐波的分布图形表⽰，这种图形称为信号的频谱图。频谱图由幅度谱图和相位谱图组成，其中幅度谱图表⽰各次谐波幅度随频率的变化关系，⽽相位谱图描述各次谐波的相位与频率的关系。根周期信号展开成傅⾥叶
级数的不同形式，频谱图⼜分为单边频谱图和双边频谱图。
  1. 单边频谱图
  根据三⾓函数形式的傅⾥叶级数展开式
f(t)=A_0+\sum_{n=1}^{\infty}A_n\cos(n\omega t+\varphi_n)
  作出的A_n与n\omega的关系称为单边幅度频谱，\varphi_n与n\omega的关系称为单边相位频谱。
  例. 画出信号f(t)=\frac{4}{\pi}[\sin(\omega_1t)+\frac{1}{3}\sin(3\omega_1t)+\frac{1}{5}\sin(5\omega_1t)+\cdots ]的单边频谱图。
  解：f(t)=\frac{4}{\pi}[\cos(\omega_1t-\frac{\pi}{2})+\frac{1}{3}\cos(3\omega_1t-\frac{\pi}{2})+\frac{1}{5}\cos(5\omega_1t-\frac{\pi}{2})+\cdots ]
  可知其基波频率为ω1，分别有⼀、三、五......次谐波分量，则其振幅谱和相位谱分别为：
  由此可见，周期函数的频谱或谱线只出现在离散点上，分布于整个频域中，形成离散谱，每条谱线间的距离为\omega_1=\frac{2\pi}{T}。离散谱是周期函数的重要特征（离散性）。此外，每⼀条谱线只能出现在基波频率ω1整数倍频率上（谐波性）。各次谐波分量的振幅虽然随着nω1的变化⽽变，但总的趋势是随着nω1的增⼤⽽逐渐减⼩（收敛性）。
  2. 双边频谱图
  根据复指数形式的傅⾥叶级数展开式
  f(t)=\sum_{n=-\infty}^{\infty}F_ne^{jn\omega t}
  设F_n=|F_n|e^{j\varphi_n}，作出的幅度|F_n|与n\omega的关系称为双边幅度频谱，相位\varphi_n与n\omega的关系称为双边相位频谱。
  例. 画出信号f(t)=\frac{4}{\pi}[\sin(\omega_1t)+\frac{1}{3}\sin(3\omega_1t)+\frac{1}{5}\sin(5\omega_1t)+\cdots ]的双边频谱图。
  根据欧拉公式：
  \sin\beta=\frac{-j}{2}(e^{j\beta}-e^{-j\beta})=\frac{1}{2}(e^{j\beta}e^{-j\frac{\pi}{2}}-e^{-j\beta}e^{j\frac{\pi}{2}})
  可得：
\begin{align*} F_n=\left\{\begin{matrix}0,\quad &n=0,\pm2,\pm4,\pm6,\cdots \\ \frac{2}{n\pi}e^{-j\frac{\pi}{2}},\quad &n=1,3,5,7,\cdots \\ \frac{2} {n\pi}e^{j\frac{\pi}{2}},\quad &n=-1,-3,-5,-7,\cdots \end{matrix}\right. \end{align*}
  则幅值为：|F_n|=\frac{2}{n\pi}\quad (n\ne0)，相位为：\varphi_n=-\frac{\pi}{2}\quad (n>0),\varphi_n=\frac{\pi}{2}\quad (n<0)
  可以看出双边谱在正负频率位置上的幅度为单边幅度谱中对应频率分量幅度的⼀半。双边谱中负频率出现仅为便于数学运算，没有物理意义，只有将负频率项与相应的正频率项合并，才是实际的频谱函数。对相位谱来说：双边相位谱以单边相位谱为基础构成奇函数。
⾮周期信号的傅⾥叶变换
  在⼯程技术中，周期函数可以展开成傅⾥叶级数，那么⾮周期函数呢?能否⽤⼀个周期函数逼近⼀个⾮周期函数呢?⼀般⽽⾔，任何⼀个⾮周期函数f(x)都可以看成是由某个周期函数f T(x)当周期T→+∞时转化⽽来的。为说明这⼀点，作周期为T的函数f T(x)使其在[-\frac{T}{2},\frac{T}{2}]之内等于
f(x)，⽽在[-\frac{T}{2},\frac{T}{2}]之外按周期T的函数f T(x)延拓出去。
f_T(x)=\left\{\begin{matrix}f(x),x \in[-\frac{T}{2},\frac{T}{2}] \\ f_T(x+T),x \notin[-\frac{T}{2},\frac{T}{2}] \end{matrix}\right.
  T越⼤，f T(x)与f(x)相等的范围也越⼤，这表明当T→+∞时，周期函数f T(x)便可以转化为f(x)。周期T趋于⽆穷⼤，其频谱间隔将趋于⽆穷⼩，从⽽信号的频谱密集成为连续频谱。同时，各频率分量的幅度也趋近于⽆穷⼩，不过这些⽆穷⼩量之间仍保持⼀定的⽐例关系。为了表明⽆穷⼩的振幅之间的差别，引⼊⼀个新的量成为“频谱密度函数”。
离散傅⾥叶变换
  从时域的采样数据转变为频域的算法，称为离散傅⾥叶变换（DFT）。DFT将采样信号的时域和频域联系起来，DFT⼴泛应⽤于谱分析、应⽤⼒学、光学、医学图像、数据分析、仪器及远程通信等⽅⾯。
  假设有N个时域采样信号，对这N个样本进⾏DFT变换，结果仍将为N个样本，但它却是频域表⽰法。时域的N个样本与频域的N个样本之间的关系如下：
  假设信号采样率为f_s，采样间隔为\Delta t，有\Delta t=1/f_s，采样信号表⽰为x_n,0\leqslant n \leqslant N-1（即有N个样本），对这N个样本进⾏傅⾥叶变化，公式如下：
X_k=\sum_{n=0}^{N-1}x_n e^{-j2\pi k n/N}\quad\quad k=0,1,2,...,N-1
  其结果X_k为x_n对应的频域表⽰。注意时域和频域中均有N个样本，同时域中的时间间隔对应的频域间隔\Delta f为：
\Delta f = \frac{f_s}{N}=\frac{1}{N\Delta t}
  \Delta f也被称为频率分辨率，增加采样次数N或减⼩采样频率f_s均能减⼩\Delta f（提⾼分辨率）
DFT计算举例
  下例为 N个采样点经DFT后的准确频率⼤⼩。假设X_0表⽰直流成分或信号的均值，为了便于观察对波形进⾏DFT后的结果，假设这个直流成分的幅度值为常数+1V，下图中列出了4个采样信号。
  每个采样点的幅值均为+1V，按时间顺序排列：x_0=x_1=x_2=x_3=1，⽤DFT公式对这个序列进⾏离散傅⾥叶变换。
  利⽤：e^{-j\theta}=\cos{\theta}-j\sin{\theta}就可以得到：
  X_0=\sum_{n=0}^{N-1}x_n e^{-j2\pi n\times 0/N}=x_0+x_1+x_2+x_3=4
  X_1=x_0+x_1(\cos{\frac{\pi}{2}}-j\sin{\frac{\pi}{2}})+x_2(\cos\pi-j\sin\pi)+x_3(\cos{\frac{3\pi}{2}}-j\sin{\frac{3\pi}{2}})=1-j-1+j=0
  X_2=x_0+x_1(\cos\pi-j\sin\pi)+x_2(\cos2\pi-j\sin2\pi)+x_3(\cos 3\pi-j\sin 3\pi)=1-1+1-1=0
  X_3=x_0+x_1(\cos{\frac{3\pi}{2}}-j\sin{\frac{3\pi}{2}})+x_2(\cos 3\pi-j\sin 3\pi)+x_3(\cos{\frac{9\pi}{2}}-j\sin{\frac{9\pi}{2}})=1+j-1-j=0
  因此，除了直流成分X_0外，其它值均为0。X_0的值取决于采样次数N，由于N=4，所以X_0=4，若N=10，则X_0=10。其它的值X_k也同样依赖于N的⼤⼩；因此为了得知频率成分的⼤⼩，经常将DFT的结果除以N。
  在Mathematica中使⽤函数可以很⽅便的计算离散傅⾥叶变换：
计算离散傅⾥叶变换
  DFT公式为:X_k=\sum_{n=0}^{N-1}x_n e^{-j2\pi k n/N}，从形式上看它是⼀个线性运算：向量\vec{x}的矩阵乘法（利⽤矩阵M将其变换到频域空间）
\vec{X}=M \vec{x}
  矩阵M(N⾏N列)可以表⽰为：
  M_{kn}=e^{-j2\pi kn/N}
  这么想的话，我们可以简单地利⽤矩阵乘法计算DFT：
import numpy as np
def DFT(x):
"""Compute the discrete Fourier Transform of the 1D array x"""
x = np.asarray(x, dtype=float)
N = x.shape[0]      # the number of samples
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x)
test_data =[1,1,1,1]
print DFT(test_data)
  结果如下：
  我们可以对⽐numpy中的FFT函数：
  可以看出FFT计算出了相同的结果，但是两者所花费的时间差别巨⼤（在ipython中进⾏测试）：
%timeit  %run  dft.py
%timeit  np.fft.fft(x)
  测试结果为：
100 loops, best of 3:  2.84ms per loop
100000 loops, best of 3:  5.09us per loop
  可以看出这种简单的DFT计算⽅法，⽐FFT慢了⼏个数量级。⼀般来说x_n和 e^{-j2\pi k n/N}都是复数，因此每计算⼀个X_k，需要N次复数乘法和N-1次复数加法，⽽k取值从0到N-1，所以完成整个DFT运算总共需要N2次复数乘法和N(N-1)次复数加法。在这些运算中乘法运算要⽐加法运算复杂，需要的时间也多⼀些。因为复数运算实际上是由实数运算来完成的，这时DFT计算公式可写为：
X_k=\sum_{n=0}^{N-1}(a+jb)(c+jd)
的形式，由此可见，1次复数乘法需要4次实数乘法和2次实数加法；1次复数加法需要2次实数加法。因⽽每运算⼀个X_k需要4N次实数乘法和
2N+2(N-1)=2(2N-1)次实数加法。所以整个DFT运算总共需要4N2次实数乘法和2N(2N-1)次实数加法。
  从上⾯统计可以看到，直接计算DFT，乘法次数和加法次数都是与N2成正⽐的，当N很⼤时，运算量是很⼤的，有时甚⾄⽆法忍受。例如对⼀幅N×N的图像进⾏DFT变化，当N=1024时，直接计算DFT所需复数乘法次数为1012次，如果⽤每秒做⼀千万次复数乘法的计算机，即使不考虑加法运算时间，也需要近28个⼩时。对实时性要求很⾼的信号处理来说，只有改进DFT计算⽅法，减少复数乘法、加法的运算次数。
快速傅⾥叶变换（FFT）
  快速傅⾥叶变换 (fast Fourier transform), 即利⽤计算机计算离散傅⾥叶变换（DFT)的⾼效、快速计算⽅法的统称，简称FFT。快速傅⾥叶变换是1965年由J.W.库利和T.W.图基提出的。采⽤这种算法能使计算机计算离散傅⾥叶变换所需要的乘法次数⼤为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。对于长度为N的输⼊⽮量，FFT是O(N logN)级的，⽽普通DFT算法是O(N^2)级的。
  FFT是怎么快速计算的呢？答案就在于它利⽤了对称性。从DFT计算公式可看出，不管输⼊信号x_n是实数还是复数，X_k总是复数(也可能虚数部分为零)，它包含两部分：幅值和相位。可以证明，对于实信号x_n，其DFT具有对称的性质：
|X_k|=|X_{N-k}|
phase(X_k)=-phase(X_{N-k})
  即幅值偶对称，相位奇对称（下标为k和N-k的两个复数共轭）。偶对称指信号关于y轴对称，奇对称指信号关于原点对称，如下图所⽰：
  下⾯的例⼦演⽰了这⼀个规律，先以rand随机产⽣有8个元素的实数数组x，然后⽤fft对其运算之后，观察其结果为8个复数：
  可以看出第(1、7)，(2、6)，(3、5)个复数互为共轭复数。DFT变换后的复数数组的另⼀规律是：下标为0和N/2的两个复数的虚数部分为0。下⾯让我们来看看FFT变换之后的那些复数都代表什么意思。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。