傅里叶变换

此篇博文主要简述傅里叶变换的相关概念以及如何代码实现离散序列的傅里叶变换。

傅里叶变换的一般形式

我们都知道，傅里叶变换是频域分析的重要工具；其将信号从时域转换到了频域，以更直观的角度向我们展示了信号的本质——频率。下面先不加证明地引出傅里叶变换的公式：

连续非周期信号
傅里叶变换为：
$$X(w)={\int }_{-\infty }^{\infty }x(t)e^{-jwt}dt$$
傅里叶反变换为：
$$x(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }X(jw)e^{jwt}dw$$

连续周期信号
傅里叶变换为：
$$X\left(n{w}_0\right)=\frac{1}{T}{\int }_{t_0}^{t_0+T}\tilde{x}(t)e^{-jn{w}_{0}t}dt$$
傅里叶反变换为：
$$x(t)=\sum _{\mathbf{n}=-\infty }^{\infty }X\left(n{w}_0\right)e^{jn{w}_{0}t}$$

其中$w_0=\frac{2\pi }{T}$

离散非周期信号
傅里叶变换为：
$$\tilde{X}\left(e^{j\Omega T}\right)=\sum _{n=-\infty }^{\infty }x(nT)e^{-jn\Omega T}$$
傅里叶反变换为：
$$x(nT)=\frac{T}{2\pi }{\int }_{-\frac{\pi }{T}}^{\frac{\pi }{T}}X\left(e^{j\Omega T}\right)e^{jn\Omega T}d\Omega$$

其中：$\Omega =k{\Omega }_0=\frac{2\pi }{NT}k=\frac{2\pi f_s}{N}k=2\pi f$，$f_s$为采样频率。

离散周期信号
傅里叶变换为：
$$\tilde{X}(k)=\sum _{n=0}^{N-1}\tilde{x}(n)e^{-j\frac{2\pi }{N}nk}$$
傅里叶反变换为：
$$\tilde{x}(n)=\frac{1}{N}\sum _{k=0}^{N-1}\tilde{X}(k)e^{j\frac{2\pi }{N}nk}$$

根据上述变换的形式，我们可以总结出如下结论：

结论

1. 连续非周期信号的傅里叶变换为非周期连续函数——傅里叶变换；
2. 连续周期信号的傅里叶变换为非周期离散函数——傅里叶级数；
3. 离散非周期信号的傅里叶变换为周期连续函数——离散时间傅里叶变换；
4. 离散周期信号的傅里叶变换为周期离散函数——离散傅里叶级数。

即时域上的离散对应频域上的周期，时域上的周期对应频域上的离散；

拓展1：对时域中的信号进行采样转化为离散序列时，相当于冲激采样信号与被采样信号时域上相乘(等价于频域上的卷积)；从而时域上的采样会引起频域上的周期延拓(这回答了为何离散化的时域信号其频谱呈现周期性)。下图详细描述了这一点：

此外从上图可以看出，由于$w_s=\frac{2\pi }{T_s}$，当$w_s<2w_m$时，采样后的信号会出现频率混叠；这要求我们要想通过采样不失真地复现原信号，采样频率必须高于信号中最高频率分量的两倍(奈奎斯特采样定理)。
拓展2：为了从采样信号中复现原信号，观察采样后的频谱图可知，只需将采样后的信号通过理想低通滤波器即可，但由于理想低通滤波器的物理不可实现性，工程中常用零阶保持器来代替；至此原信号通过冲激采样信号采样再通过零阶保持器便完成了信号的采样和保持，最终得到了阶梯状的信号序列，如下图所示：

离散傅里叶级数

接着来看离散周期信号的傅里叶变换(即离散傅里叶级数$DFS$)
其傅里叶变换如下所示：
$$\tilde{X}(k)=\sum _{n=0}^{N-1}\tilde{x}(n)e^{-j\frac{2\pi }{N}nk}$$
实际上，由于离散信号的傅里叶变换是周期函数，故$\tilde{X}(k)$有无穷多项，为方便计算，取$k$的上限为$N-1$(事实上，此时系数矩阵满秩，方程具有唯一解)；将上式写成矩阵形式则有：
$$\tilde{X}=W_N\tilde{x}$$
即：
$$\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& W_N^1& \cdots & W_N^{N-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& W_N^{N-1}& \cdots & W_N^{(N-1{)}^2}\end{array}\right]\left[\begin{array}{c}\tilde{x}(0)\\ \tilde{x}(1)\\ ⋮\\ \tilde{x}(\mathrm{N}-1)\end{array}\right]=\left[\begin{array}{c}\tilde{X}(0)\\ \tilde{X}(1)\\ ⋮\\ \tilde{X}(\mathrm{N}-1)\end{array}\right]$$

其中：

$$W_N=e^{-j\frac{2\pi }{N}}$$

上述矩阵即为DFS矩阵，观察其形式可知其为范德蒙矩阵，则其逆可表示为：

$$W_N^{-1}=\frac{1}{N}{W}_N^{\ast }$$

进一步可以写成：
$$W_N^{-1}={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& W_N^1& \cdots & W_N^{N-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& W_N^{N-1}& \cdots & W_N^{(N-1{)}^2}\end{array}\right]}^{-1}=\frac{1}{N}\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& W_N^{-1}& \cdots & W_N^{-(N-1)}\\ ⋮& ⋮& \ddots & ⋮\\ 1& W_N^{-(N-1)}& \cdots & W_N^{-(N-1{)}^2}\end{array}\right]$$
上述步骤的详细证明

同样地，可以写出离散周期信号的傅里叶反变换形式为：

$$\tilde{x}=W_N^{-1}\tilde{X}$$

傅里叶变换在MATLAB中的实现

上述步骤在Matlab中的实现为：

%DFS：
function [Xk] = dfs(xn,N)
n=0:1:N-1;
k=0:1:N-1;
WN = exp(-j*2*pi/N);
nk = n'*k;
WNnk = WN.^nk;
Xk = xn*WNnk;
end

%IDFS:
function [Xk] = dfs(xn,N)
n=0:1:N-1;
k=0:1:N-1;
WN = exp(-j*2*pi/N);
nk = n'*k;
WNnk = WN.^(-nk);
Xk = xn*WNnk/N;
end

一个例子：

Ts = 0.01; %离散步长
T = 10;    %采样时间
t= Ts:Ts:T;%采样步长
w1 = 5*2*pi;
w2 = 20*2*pi;
w3 = 35*2*pi; %频率分量，单位为弧度
N=length(t);  %样点个数
y= cos(w1*t) + cos(w2*t) + cos(w3*t);
Fs=1/Ts;                 %采样频率
df=Fs/(N-1) ;            %分辨率
f=(0:N-1)*df;            %其中每点的频率
Xf = dfs(y,N);           %傅里叶变换
Phase = zeros(1,len);    %用于存储各频率点的幅值
Amplitude = zeros(1,len);%用于存储各频率点的相位
%幅值
for i=1:N
    Amplitude(i) = abs(Xf(i));
    Phase(i) = angle(Xf(i));
end
%绘制源信号波形
figure(1)
plot(t,y)
grid on
%绘制频谱图
figure(2)
plot(f,2*Amplitude/N)
xlabel('频率/Hz')
ylabel('幅值')
grid on

原信号波形：

频谱图(原图)：

频谱图(截断处理后):

注意
1.对于一个时域信号$x(t)$，采样频率为$F_s$，采样点数为$N$，采样后序列为$y(n)$,则上述$F=dfs(y,N)$计算所得的结果就是一个有$N$个元素的复数序列$a_n+b_{n}j$；该序列每一个复数对应着一个频率点，复数的模值和相角对应信号在该频率值下的幅相特性；
2.频率：$F$序列的第$n$点处的对应的实际频率是$n\ast Fs/N$；
3.幅值：$F$序列的第$n$点处对应的复数为$a+bj$，模值$A=\sqrt{a^2+b^2}$，那么实际信号在该频率值下的幅值是$2\ast A/N$;特别的，当$n=0$时，对应幅值为$A/N$；
4.相位：$F$序列每个频率点处的复数的相位就是信号在该频率值下的相位；
5.采样频率必须大于信号中最高频率分量的两倍，否则会引起频率混叠；
6.$DFS$得到的$F$序列关于中间位置$F_s/2$对称(周期性)，由采样定理可知，高于$F_s/2$的信号将会出现频率混叠，频谱图显示的最高频率应为$F_s/2$，故需对$F$序列的横坐标进行截断处理。

离散傅里叶变换

离散傅里叶级数(DFS)提供了一种对离散时间傅里叶变换做数值计算的技巧，它在时域和频域都是周期的，但在实际中大多数信号不具有周期性，且很可能只具有有限持续时间。这就要求我们探讨一种对非周期离散有限信号序列进行傅里叶变换的方法，即离散傅里叶变换(DFT)。
DFT是对任意有限持续时间序列可数值计算的傅里叶变换
长度为N的有限长序列$x(n)$可以看作周期序列$\tilde{x}(n)$的主值序列，同样的$\tilde{x}(n)$可以看作$x(n)$的周期延拓。

由于无频率混叠，则时域中的一个周期长的主值序列对应于频域中的一个周期长的主值序列。从DFS的时域和频域中各取出一个周期，即可得到有限长度的离散序列的时域和频域傅里叶变换。

$$\{\begin{array}{l}X(k)=\left[{\sum }_{n=0}^{N-1}\tilde{x}(n)W_N^{nk}\right]R_N(k)=\tilde{X}(k)R_N(k)\\ x(n)=\left[\frac{1}{N}{\sum }_{k=0}^{N-1}\tilde{X}(k)W_N^{-nk}\right]R_N(n)=\tilde{x}(n)R_N(n)\end{array}$$
其中$R_N(k)$和$R_N(n)$均为窗函数，负责从时域与频域序列中取出主值序列。即可得到：

$$\{\begin{array}{l}X(k)=DFT[x(n)]={\sum }_{n=0}^{N-1}x(n)W_N^{nk}0\le k\le N-1\\ x(n)=\mathrm{IDFT}[X(k)]=\frac{1}{N}{\sum }_{k=0}^{N-1}X(k)W_N^{-nk}0\le n\le N-1\end{array}$$

事实上，DFS与DFT的表达式没有本质区别。