[ MATLAB ] 傅里叶变换（三）：傅里叶变换

前言

专题的前两篇文章（ [ MATLAB ] 傅里叶变换（二）：傅里叶级数（复指数表示）），我们讨论了连续周期信号傅里叶级数的两种表示形式，初步建立了频谱的概念。然而，就实际经验而言，非周期信号才是主流。因此，这篇文章将讨论非周期连续信号的谱密度（通常简称为频谱），即大名鼎鼎的傅里叶变换FT，并用Matlab仿真加强理解。

理论回顾

首先，我们借助傅里叶级数$\to$傅里叶变换的过程，回顾一下傅里叶变换的基本数学知识。

若$x(t)$是周期连续信号，周期为$T$，模拟角频率${\Omega }_0=\frac{2\pi }{T}$，则$x(t)$的傅里叶级数为：
$$正变换:X\left(jk{\Omega }_0\right)=\frac{1}{T}\underset{0}{\overset{T}{\int }}x\left(t\right)e^{-jk{\Omega }_{0}t}dt$$
$$反变换:x\left(t\right)=\sum _{k=-\infty }^{\infty }X\left(jk{\Omega }_0\right)e^{jk{\Omega }_{0}t}$$
现在考虑$x(t)$为非周期信号，即$T\to \infty$，此时${\Omega }_0\to 0$，${\Omega }_0$可视作相邻谱线的频率间隔$(k+1){\Omega }_o-k{\Omega }_0$，因此，$k{\Omega }_o=\Omega$，代入正变换可得
$$X{\left(j\Omega \right)}^{\prime }=\frac{1}{T}\underset{-\infty }{\overset{\infty }{\int }}x\left(t\right)e^{-j\Omega t}dt$$
显然，$X{\left(j\Omega \right)}^{\prime }$是一个无穷小量，即$x(t)$的频率分量振幅无穷小。但是，无穷小之间仍然有相对大小。令
$$X\left(j\Omega \right)=TX{\left(j\Omega \right)}^{\prime }=\underset{-\infty }{\overset{\infty }{\int }}x\left(t\right)e^{-j\Omega t}dt$$
$X\left(j\Omega \right)$即为频率分量的相对大小，称其为$x(t)$的频谱密度（通常简称为频谱）。上式同样被称为FT的正变换，给出完整的FT正反变换表达式：
$$正变换:X\left(j\Omega \right)=\underset{-\infty }{\overset{\infty }{\int }}x\left(t\right)e^{-j\Omega t}dt$$

$$反变换:x\left(t\right)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}X\left(j\Omega \right)e^{j\Omega t}d\Omega$$
这一小节，从数学角度简单分析了FT的由来。相信很多铁子们仍然对谱密度这个概念有些困惑，这将在下面的Matlab仿真中进行梳理验证。

Matlab仿真常见信号的频谱

这里算是开胃小菜，我们用Matlab仿真两个非周期信号的傅里叶变换。

（1）连续时间信号$x\left(t\right)=t{e}^{-|t|}$，求出其频谱并绘制幅度谱；

（2）某信号的频谱表达式为$X\left(j\Omega \right)=\pi e^{-|\Omega |}$ ，求出其时域表达式并绘制图像；

syms t w;  

f1 = t*exp(-abs(t));  % 时域表示
subplot(2,2,1);ezplot(f1);
F1 = fourier(f1);  % 正变换
subplot(2,2,2);ezplot(abs(F1));

F2 = pi*exp(-abs(w));  % 频域表示
subplot(2,2,3);ezplot(abs(F2));
f2 = ifourier(F2);  % 反变换
subplot(2,2,4);ezplot(f2);

Matlab探索谱密度概念

前面的公式推导可以看出，单个频率分量对应的幅值趋近于0，并且频率分量数量趋近于无穷，因此我们用傅里叶变换得到的是谱密度。

谱密度完全可以类比于物理中的密度。举例说明，想象有一块石头，我们很容易算出他的密度是多少，但是我们取石头上非常非常小的一块，它的质量趋于0，它的体积也趋于0。这时就要引出密度的概念了，无穷小的质量/无穷小的体积得到的常数，就是这一小块的密度了。

现在，我们想知道石头中某一小块的质量，即便这一小块体积非常非常小，还是需要测量某一体积的石头对应的质量。因此，如果需要获得某一频率分量对应的幅值，必须要用该分量的谱密度值乘以采样频率区间宽度。采样频率区间宽度类比于某一小块石头的体积，分量幅值类比于质量，谱密度类比于密度。
下面的仿真中，假设矩形脉冲信号的谱密度已知；用分量的谱密度值乘以采样频率区间宽度，得到分量的幅值；最后叠加所有分量，还原矩形脉冲信号。

• step 1 划分采样区间和采样间隔

syms t;  
f = rectangularPulse(t);  % 矩形脉冲信号，原信号

st = 0; ed = 500;  % 频域采样区间
components_num = 5000;  % 频域采样项数
freqs = zeros(1,components_num);  % 分量频率/采样频率
width = (ed-st)/components_num;  % 区间宽度
sample_intervals = linspace(st,ed,components_num+1);  % 分割采样区间
for idx=1:length(sample_intervals)-1
    freq = sample_intervals(idx)+(sample_intervals(idx+1)-sample_intervals(idx))/2;
    freqs(idx) = freq;                   
end

step 2 计算频谱值并乘以频率区间，从而得出分量幅值，叠加近似

F = fourier(f);  % 傅里叶变换
F_function = matlabFunction(F);
components = sym('c',[1,components_num]);  % 分量
idx = 1;
for freq=freqs
    amp_weight = F_function(freq);  % 幅值密度
    amp = 2*abs(amp_weight)*width/2/pi;  % 密度+搬移, 注意除以2pi（反变换表达式）
    ang = angle(amp_weight);
    component = amp*cos(freq*t+ang);  
    components(idx) = component;
    idx = idx+1;
end

随着我们采样的频率范围越大，采样数量（分量数量）越多，叠加得到的信号愈发逼近原始的矩形脉冲信号，验证成功。从上面仿真中，我们也可以发现，谱密度也可以理解为权重，即谱密度越大，乘以相同采样间隔后，得到的分量幅值也就越大。

仿真总结

1. 可以采用物理中的密度的方式类比谱密度的概念，从而理解傅里叶变换中谱密度的意义。
2. 不需要再执着于分量幅值的绝对大小，而是聚焦于相对大小。
3. 目前我们成功的实现了连续非周期时间信号的频谱变换，但是计算机运算中，仍然喜欢离散的有限值，下面我们将对连续时间信号进行离散化，引出离散时间傅里叶变换的概念（DTFT）的概念，由与频域和时域均是离散的目标进了一步。

创作不易，希望铁子们可以给出宝贵的建议，这将是我这个小菜猫不断前进的动力，谢谢！