DFT

对DFT相关的内容做一个梳理。

傅里叶变换（FT）

在学习信号处理的过程中，我们最早接触的是时域连续信号的傅里叶变换：
$X(j\Omega) = \int_{-\infty}^{+\infty} x(t)e^{-j\Omega t}dt$X(jΩ)=∫−∞+∞​x(t)e−jΩtdt

而在实际应用中，经常接触到的是序列，即时域离散信号。我们可以把上式中的积分变为求和，自然而然地得到时域离散信号的傅里叶变换：
$X(e^{j\omega}) =\sum_{n = -\infty}^{+\infty}x(n)e^{-j\omega n}$X(ejω)=n=−∞∑+∞​x(n)e−jωn

需要注意的是，$X(e^{j\omega})$X(ejω)具有$2\pi$2π周期性。$\omega=0$ω=0处代表低频，$\omega=\pi$ω=π处代表高频（这个下面会有解释）。

离散傅里叶变换（DFT）

上述时域离散信号的FT在频域上是连续的，这样对于计算机来说还是不能处理。这时候，DFT诞生了。设序列$x(n)$x(n)长度为M，它的N点DFT为：
$X(k) =\sum_{n = 0}^{N - 1}x(n)e^{-j\frac{2\pi}{N}kn},\quad k=0, 1, ..., N-1$X(k)=n=0∑N−1​x(n)e−jN2π​kn,k=0,1,...,N−1

要求$N\ge M$N≥M（原因见下文）。可以看出，它刚好是有限长序列$x(n)$x(n)的傅里叶变换$X(e^{j\omega})$X(ejω)在$[0, 2\pi]$[0,2π]上的等间隔采样。（因为FT具有$2\pi$2π周期性，所以在$[0, 2\pi]$[0,2π]上采样就够了）

采样定理

时域采样定理

时域采样定理是这样的：

对模拟信号$x_a(t)$xa​(t)进行时域等间隔理想采样（冲激采样），则理想采样信号$\hat x_a(t)$x^a​(t)的频谱是原模拟信号频谱以采样频率为周期的周期延拓。对于带限信号，当采样频率大于2倍信号频率时，可以由采样信号无失真恢复原模拟信号。

上述采样定理针对的是理想采样，这在现实中不可实现。实际中采用的是时域离散采样，即$x(n) = x_a(nT)$x(n)=xa​(nT)（注意：理想采样$\hat x_a(t)$x^a​(t)仍是模拟信号，而$x(n)$x(n)是时域离散信号）。$x(n)$x(n)的傅里叶变换$X(e^{j\omega})$X(ejω)与$x_a(t)$xa​(t)的傅里叶变换$X(j\Omega)$X(jΩ)的关系为：
$X(e^{j\omega})=1/T\times\sum_{k=-\infty}^{+\infty}X_a(j\Omega-jk\Omega_s)$X(ejω)=1/T×k=−∞∑+∞​Xa​(jΩ−jkΩs​)

其中，$\Omega_s=2\pi/T=2\pi F_s$Ωs​=2π/T=2πFs​，$\omega=\Omega T$ω=ΩT。

$X(e^{j\omega})$X(ejω)是$\omega$ω的函数，$X_a({j\Omega})$Xa​(jΩ)是$\Omega$Ω的函数，二者之间通过一个坐标变换 $\omega=\Omega T$ω=ΩT 联系在了一起。也就是说，时域离散采样信号的频谱是模拟信号频谱以$\Omega_s$Ωs​为周期做周期延拓之后，再在横轴上做一个伸缩变换，再乘以一个系数得到的。

同理，对于时域离散采样，也要满足采样频率大于2倍信号频率，才能避免频谱混叠。

根据$\omega=\Omega T$ω=ΩT，$\omega=2\pi$ω=2π对应于$\Omega=\Omega_s$Ω=Ωs​，$\omega=\pi$ω=π对应于$\Omega=\Omega_s/2$Ω=Ωs​/2。$\Omega_s/2$Ωs​/2是折叠频率，这就是为什么说$\omega=\pi$ω=π对应于高频部分了。

进一步地，对$x(n)$x(n)做DFT，就相当于对原模拟信号频谱的周期延拓在$[0, \Omega_s]$[0,Ωs​]上等间隔采样！

频域采样定理

对任意序列$x(n)$x(n)的傅里叶变换$X(e^{j\omega})$X(ejω)，在$[0, 2\pi]$[0,2π]上等间隔采样，得到$\tilde X_N(k)$X~N​(k)（这里不限制k的取值范围）。可以知道，$\tilde X _N(k)$X~N​(k)是以N为周期的序列。根据离散傅里叶级数理论，$\tilde X_N(k)$X~N​(k)必然是某个周期序列$\tilde x_N(n)$x~N​(n)的DFS。通过推导可知，$\tilde x_N(n)$x~N​(n)是原序列$x(n)$x(n)以N为周期的周期延拓。分别取$\tilde X_N(k)$X~N​(k)与$\tilde x_N(n)$x~N​(n)的主值序列$x_N(k)$xN​(k)和$x_N(n)$xN​(n)，二者构成一对DFT。

从而可以得到频率采样定理：设序列$x(n)$x(n)长度为M，对$X(e^{j\omega})$X(ejω)在$[0, 2\pi]$[0,2π]上N点等间隔采样（也就是做N点DFT），仅当 N>=M 时，才能由$X_N(k)$XN​(k)恢复$X(n)$X(n)，避免混叠。

对比：
时域采样定理要求信号频率有限；频率采样定理要求序列长度有限。

MATLAB中fft的使用

以下代码摘自MATLAB help文档

一、定义信号，注意几个关键量：

Fs = 1000;            % 采样频率   
T = 1/Fs;             % 采样间隔
L = 1500;             % 序列长度
t = (0:L-1)*T;        % 时间向量
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 		% 两个正弦信号
X = S + 2*randn(size(t));			% 加噪声

二、傅里叶变换

Y = fft(X);			% 得到的是原信号频谱在[0, Fs]上的L点等间隔采样
P2 = abs(Y/L);				% 这个为什么要除以L，请看下面解释
P1 = P2(1:L/2+1);		% 实信号的FT是共轭对称的，我们只看半边！
P1(2:end-1) = 2*P1(2:end-1);	% 由于是单边谱，所以幅度乘2
f = Fs*(0:(L/2))/L;				% 频率轴，采样间隔是Fs/L，频率范围是[0, Fs/2]
plot(f,P1) 

这里有一个疑问是：为什么要对fft得到的结果除以L，按照时域采样定理那一部分的结论：
$X(e^{j\omega})=1/T\times\sum_{k=-\infty}^{+\infty}X_a(j\Omega-jk\Omega_s)$X(ejω)=1/T×k=−∞∑+∞​Xa​(jΩ−jkΩs​)

我们似乎应该对fft的结果乘以T，也就是除以Fs，才能得到原频谱。找了谷歌，明白了原因：

根据Parseval定理，信号的能量对于其频谱的能量，即：
$\int_{-\infty}^\infty | x(t) |^2 \, dt = \int_{-\infty}^\infty | X(f) |^2 \, df$∫−∞∞​∣x(t)∣2dt=∫−∞∞​∣X(f)∣2df

DFT是对频谱的L点采样：
$X_k = Fs\times X(\frac{F_s}{L}\times k)，k =0,...,L-1$Xk​=Fs×X(LFs​​×k)，k=0,...,L−1

我们用求和来近似积分：

$\int_{-\infty}^\infty | X(f) |^2 \, df\approx \sum_{k=0}^{L-1}X(\frac{F_s}{L}\times k)\Delta f =\frac{1}{L}\sum_{k=0}^{L-1}X_k$∫−∞∞​∣X(f)∣2df≈k=0∑L−1​X(LFs​​×k)Δf=L1​k=0∑L−1​Xk​

即：
$\sum_{k=0}^{L-1}X_k= {L}\int_{-\infty}^\infty | X(f) |^2 \, df$k=0∑L−1​Xk​=L∫−∞∞​∣X(f)∣2df

可见，采样点数越多，右边的和越大，即信号的能量被放大越多。我们希望能够在频谱上反应信号的能量，因此要对fft的结果除以L。

其实，如果从频谱大小来看，确实应该除以Fs。但我们从能量的角度来考虑，却是需要除以L！