彻底搞懂--离散傅里叶变换和频谱分析

彻底搞懂–离散傅里叶变换和频谱分析

举头望明月 低头傅里叶

1、名词解释

	时域	频域
FS（Fourier Series）	连续、周期	离散、非周期
FT（Fourier Transform）	连续、非周期	连续、非周期
DFS（Discrete Fourier Series）	离散、周期	离散、周期
DTFT（Discrete-time Fourier Transform）	离散、非周期	连续、周期
DFT（Discrete Fourier Transform）	离散、周期	离散、周期
FFT（Fast Fourier Transform）	离散、周期	离散、周期

对于连续时间的傅里叶变换，不再进行推导和说明，工程中应用到的主要是DFT和FFT，因为计算机只能处理离散的数据，为此对离散时间的傅里叶变换做重点说明：

• DFS 可以处理离散周期时间序列，但是周期性数据是无穷大的，虽然时间是离散的，但是有无穷多个周期，计算机无法处理无穷多的数据，对于频域的数据也是一样。
• DTFT 可以处理离散非周期时间序列，虽解决了时域问题，但变换到频域之后，频域仍是连续的，计算机无法处理连续的频域数据。
• DFT DFT处理的是离散、周期的时间序列，但是这个“周期”需要加一个引号，因为这是人为造的周期（即周期延拓），而非原始信号的周期，频域上也是如此。这样计算机就可以在无穷多个周期中，任取一个周期进行处理。
• 总结：原始信号$\Rightarrow$时域采样$\Rightarrow$时域截断$\Rightarrow$频域采样，得到离散的、有限的、计算机可以处理的数据。

2、离散傅里叶变换（DFT）

1. 原始信号进行傅里叶变换：

2. 原始信号周期性采样后的傅里叶变换：时域周期采样，频域周期延拓

3. 频域周期性采样后的反傅里叶变换：频域周期采样，时域周期延拓

4. 原始信号周期性采样、频域周期性采样，综合（b）和（c）：时域周期延拓、频域周期延拓

以上四个步骤的数学表达为：

$f(t)$乘上一个周期性冲激序列得到$f_p(t)$，表现为离散化采样，$\stackrel{-}{f_p(t)}$表示$f_p(t)$与一个周期性冲激序列的卷积，表现为周期延拓。对应的频域操作与之类似。

需要注意的是FT的性质：时域上卷积，对应于频域上相乘；而频域上的卷积等价于时域上的相乘（差一个参数1/2$\pi$），另外还要注意采样频率，要满足奈奎斯特定理，即采样频率$f_s\ge 2f_{max}$(原始信号的最大频率)，一般要大于2倍。

3、DFT推导

1. $f_p(t)\Rightarrow F_p(jw)$

2. ${\stackrel{-}{f}}_p(t)\Rightarrow \stackrel{-}{F_p(jw)}$：注意下面用到了FT的卷积性质，即时域卷积，频域相乘；$F_p(jw)$和$\Omega {\sum }_{n=-\infty }^{\infty }\delta (\omega -m\Omega )$分别为$f_p(t)$和周期性冲激序列的傅立叶变换，代入上面$F_p(jw)$的值：

其中：
$\Omega =2\pi /T$：T为采样所用的总时间，有$T=N\ast T_s$ ，即在T时间段内，以$T_S$间隔，共采样N个点，则$\Omega T_s=2\pi /N$
m为（1，2，3，…）对频域离散采样，得到的索引（序号）根据冲激信号的性质，只有当$\omega =m\Omega$时才有非零值，其余点均为零（如上图d所示）
3. 代入$\Omega T_s=2\pi /N$

同时与上图9-1中的：

对比化简得到：

将频域上的索引m，和时域上的索引k看作自变量；频域上以${\omega }_s$的频率大小为周期，只需在一个周期内取N个点（截取），就可以表示全部（如上图d所示）

得到最终的DFT（离散傅里叶变换）表达式：$W_N^{mk}=e^{-j\frac{2mk\pi }{N}}$

用相似的方法，根据IFT的性质，用上述方法可以推导出IDFT（离散傅里叶反变换）表达式：

4、Matlab频谱分析

T=10.00;%总时间
f=1/T;%f可理解为频域上的采样频率间隔
Ts=0.01;%采样时间间隔
fs=1/Ts;%时域上的采样频率
t=0.01:Ts:10;%表示从0.01到10，间隔为Ts，（0.00，0.01，0.02......）采样频率为100Hz,采样1000个点
signal=5*sin (2*pi*10*t+.75*pi);%一个频率为10Hz的正弦信号
mix = 15*randn(1,length(signal));%根据原信号生成一个噪声
mix_signal = signal + mix;%把原始信号复杂化，加上噪声

%傅里叶变换
fft_signal=fft(signal);%调用fft函数进行离散傅立叶变换
fft_mix_signal=fft(mix_signal);%同上
N=length(fft_signal);%得到离散点的个数

%fft后，得到的并非原始频率和幅值，需进行恢复
f=fs/N;%恢复为原始频率
x_1=(0:N-1)*f;%相当于频域上离散的采样点
x_2=(0:N-1)*f;%同上

figure(1);
plot(x_1,abs(fft_signal));%abs()对复数进行取模，画出图形
figure(2);
plot(x_2,abs(fft_mix_signal));%同上

%上述恢复了频率，但得到的并非严格的双边频谱，且未恢复幅值，需要进行移频，
shiftfft_signal=fftshift(fft_signal);%fftshift()进行移频
shiftfft_mix_signal=fftshift(fft_mix_signal);%同上
shift_x_1=(-N/2:N/2-1)*f;%此时的横坐标关于原点对称，因此以0为中心取N个点
shift_x_2=(-N/2:N/2-1)*f;%同上

figure(3);
plot(shift_x_1,abs(2*shiftfft_signal/N));%恢复幅值
figure(4);
plot(shift_x_2,abs(2*shiftfft_mix_signal/N));%同上

运行结果如下：

5、频率和幅值的复原

1. 共轭性质
   考虑$F(m)和F(N-m)$；
   $F(m)={\sum }_{k=0}^{N-1}f(k)W_N^{mk}$，其中$W_N^{mk}=e^{-j\frac{2mk\pi }{N}}$；
   $F(N-m)={\sum }_{k=0}^{N-1}f(k)W_N^{(N-m)k}$，其中$W_N^{(N-m)k}=e^{-j2\pi }{e}^{j\frac{2mk\pi }{N}}$，其中$e^{-j2\pi }=1$，通过欧拉公式推导；
   $F(N-m)={\sum }_{k=0}^{N-1}f(k)W_N^{-mk}$，所以$F(m)和F(N-m)$二者共轭；

2. 幅值
   $Nf(k)={\sum }_{m=0}^{N-1}F(m)W_N^{-mk}$；
   $Nf(k)=F(0)+F(1)W_N^{-k}+F(2)W_N^{-2k}+...+F(c)W_N^{-ck}+...+F(N-1)W_N^{-(N-1)k}$，c为中间某项；
   复数可表示为：$F(c)=a_c+j{b}_c$，$F(N-c)=a_c-j{b}_c$；

   欧拉公式：$e^{jx}=cosx+jsinx$，代入$W_N^{(N-m)k}$

   $F(c)W_N^{-ck}+F(N-c)W_N^{-(N-c)k}$
   $=2a_{c}cos\frac{2ck\pi }{N}-2a_{c}sin\frac{2ck\pi }{N}$
   $=2\sqrt{a_c^2+b_c^2}cos(\frac{2ck\pi }{N}+{\phi }_c)$

   即如下代码所示

   plot(shift_x_1,abs(2*shiftfft_signal/N));%取模、扩大2倍
   

3. 频率
   $f(k)\Rightarrow f(k{T}_S)\Rightarrow f(k/f_S)$，k表示第k个采样点，对应的时间是$k{T}_s$，最终只需把k替换成$f_{s}t$即可得到初始信号$f(t)$
   $f(k)\Rightarrow f(t)$
   根据上述幅值的推导：
   $F(c)W_N^{-ck}+F(N-c)W_N^{-(N-c)k}$
   $=2a_{c}cos\frac{2ck\pi }{N}-2a_{c}sin\frac{2ck\pi }{N}$
   $=2\sqrt{a_c^2+b_c^2}cos(\frac{2ck\pi }{N}+{\phi }_c)$
   即N项和式中，第c项、倒数第c项之和，可以用上式表达
   那么：
   $Nf(k)={\sum }_{c=0}^{N/2-1}2\sqrt{a_c^2+b_c^2}cos(\frac{2ck\pi }{N}+{\phi }_c)$
   $f(k/f_s)\Rightarrow f(t)$，把k替换成$f_{s}t$
   $Nf(t)={\sum }_{c=0}^{N/2-1}2\sqrt{a_c^2+b_c^2}cos(\frac{c{f}_s}{N}(2\pi t)+{\phi }_c)$
   从而得到$f=\frac{c{f}_s}{N}$，c取（0，1，2…N/2-1）
   即如下代码所示：

   f=fs/N;%恢复为原始频率
   x_1=(0:N-1)*f;%相当于频域上离散的采样点
   x_2=(0:N-1)*f;%同上
   

   这里对每个点都进行了恢复；
   事实上只需要一半，因为当c>N/2时，$F(m)和F(N-m)$就开始和前面(0,N/2)重复。

6、参考文献

1. 《信号与线性系统》,第五版，孟桥等著

2. 知乎-深入理解离散傅里叶变换(DFT)

7、说明

本人初出茅庐，刚刚接触信号系统，信号小白，对于DFT和频谱分析，卡了很长时间，搞懂了离散傅里叶变换，但又不懂怎么频谱分析，经过了漫长的黑夜。
本文的推导以及演示过程，算是呕心沥血，但皆为本人个人理解，若有错误和不足，欢迎指出，我们一起探讨，共同成长。