Math

从应用的角度说说傅立叶变换

Posted by MaZhaoxin on January 9, 2020

傅立叶变换给我们了一种新的视角来看待世界，今天从应用的角度简单总结一下其使用上需要注意的地方。

DFT与FFT

我们都知道，对于周期性的信号，其傅立叶变换结果为离散值，称为傅立叶级数，即任何周期信号都可以被看做不同系数的正余弦信号相加的结果。傅立叶变换的结果表征了信号在频域的特征，因此说傅立叶变换把时域信号转换为了频域信号。如果只看频域信号的幅度信息，即为信号的频谱（spectrum）。通过观察频谱，可以很方便的判断信号中各频率成分的大小。

傅立叶级数中存在关系：信号的周期与前面所述的正余弦信号的最低频率（基频）呈反比。当信号的周期趋于无穷大时，基频趋于零。由于可以把非周期信号看作周期为无穷大的信号，因此非周期信号的傅立叶变换结果为连续值。

而在实际应用中，需要通过采样来获取信号的数值，因此离散时间域傅立叶变换（DFT）更具有实用价值。至于计算方法，将连续时间域的傅立叶变换中的积分换成累加即可得到离散时间域的傅立叶变换公式。另外由于采样时所用的冲击函数会使频谱出现周期延拓，需要注意一般情况下采样频率不能低于被采样信号截止频率的一半（奈奎斯特采样定理）。

由于DFT的计算过程中需要很多次计算（对于长度为$N$的复数信号，每个采样点要进行$N$次复数乘法和$N-1$次复数加法），基2的快速傅立叶变换（FFT）算法就被发明了出来。FFT的核心思想是分段（一分为二）计算，组合输出。这种算法可以在点数较多时大大减少计算量，但要求采样点数必须是2的整数次幂。

基本性质

  1. $N$个点的输入就有$N$个点的输出；
  2. 对于采样频率为$f_s$的输入信号，输出对应的频率为$0 \sim \frac{N-1}{N}\cdot f_s $，通过频域的延拓可以转换到$-\frac{N}{2}\cdot f_s \sim \frac{N-1}{2N}\cdot f_s$；
  3. 由1、2不难看出，频谱的分辨率带宽（RBW）是$f_s/N$；
  4. 当输入信号为实数信号时，频谱关于DC对称，因此对于一般的电信号只看正频率部分的频谱即可；
  5. 对于一段固定长度为$t_{signal}$的波形，如果对其采$N$个点，可得
    • 采样周期为$T_s=t_{signal}/N$，对应的采样频率为$f_s=1/T_s=N/t_{signal}$；
    • 可以观察得到的最大频率成分在$\frac{N-1}{2N}\cdot f_s = \frac{N-1}{2}\cdot \frac{1}{t_{signal}}$，即采样点数越多，可以观察的频谱范围越广；
    • 频谱的精度用分辨率带宽（RBW）描述，有$RBW=1/t_{signal}$，即信号时间越长，频谱越精细；
    • 对于数字示波器而言，采样深度（采样点数）是定值，需要在采样信号时间长度和采样率之间做折中，即在频谱范围和频谱精度之间做折中；
  6. 由FFT直接得到的频谱表征的为每RBW频率范围内的信号幅度大小，即功率谱，将其再除以RBW可以得到功率谱密度（PSD）；
    • 当单音信号占主导时，减小RBW可以看到PSD的峰值会逐步增加；
    • 当噪声或宽带信号占主导时，减少RBW并不会影响PSD；
    • 因此可以通过不断减小RBW来判断噪底（noise floor）处于什么水平，但这是以增长采样信号长度为代价的；

频谱泄漏

当被采样信号中有明显的周期信号，且做FFT的部分没有包含完整的周期时，由于其首尾处的相位不连续导致信号频谱中各谱线会相互影响，最终使结果失真。

如果被采样信号比较简单，可以观察到明显的周期，那么选取合适的采样范围是避免频谱泄漏的最佳方案。

否则需要通过加窗函数的方式减弱首尾处的幅度。而时域上的相乘对应了频域上的卷积，因此加窗函数会使频谱出现旁瓣，在窗函数的选择上需要注意，常用的窗函数有Hanning、Hamming等。另外，在计算功率时注意归一化问题。

信号混叠

由于频谱的周期延拓，当采样频率不是被采样信号频率的整倍数时，会在频谱上出现一些混叠信号。当观察有高功率谐波的周期性信号且采样频率不够高时尤其明显，会影响对频谱纯净度的判断。

解决办法是增加采样频率（即在频谱精度和频谱范围之间折中）或选择合适的采样频率。

例如，当需要观察一个频率为$f_{clk}$的时钟信号的频谱时，如果既不想看到频谱泄漏又不想看到信号混叠，需要如何选择采样参数？

  • 时钟周期$T_{clk}$，信号最大长度$t_{max}$，做FFT的部分信号长度$t_{signal}$，采样周期$T_s$，采样频率$f_s$，采样点数$N$；
  • 为了不出现频谱泄漏，要求$K_1 = t_{signal}/T_{clk}$为整数；
  • 为了不出现信号混叠，要求$K_2 = f_s/f_{clk} = T_{clk}/t_s = N*T_{clk}/t_{signal} = N/K_1$也要为整数；
  • 由于FFT要求$N$为2的整数次幂，即$N$的素因子只有2，那么$K_1$与$K_2$也均为为2的整数次幂；
  • 所以，$t_{signal} = T_{clk} \cdot 2^{floor(\log_{2}{(t_{max}/T_{clk}}))}$，而$N$取尽量的大即可。

频谱降噪

一般情况下，FFT的结果充满了由于各种因素导致的噪声，表现在图表上即为波动范围很大的谱线，为了得到更加准确（平滑）的结果还需要对其进行降噪。

比较常用的是Welch方法（MATLAB中有pwelch函数可以调用，输出的结果是PSD）。对于总长度为$N$的信号，每次取$N_{FFT}$个点为一段，加窗函数后做FFT，相邻两段之间重复取$N_{overlap}$个点（一般取$N_{FFT}$的33~50%，当然也可以为0），最后对每段的结果求平均作为最终结果。

下面举例说明MATLAB中pwelch函数的使用方法：

function showPSD(w, fs)
% by:MaZhaoxin @2015/6/17
if nargin<2
    fs = 1;
end

% 截取w，以保证其长度满足FFT的需求
L = 2^nextpow2(length(w));
if L~=length(w)
    fprintf('\t***Warning: Samples will be cut(%d -> %d).\n', length(w), L/2);
    L = 2^(nextpow2(length(w))-1);
    w = w(end-L+1:end);
end

% 共计分为7段
NFFT = L/4;
[p, f] = pwelch(w, hann(NFFT), NFFT/2, NFFT, fs, 'twosided');

% 处理结果，将频率映射到-fs/2 ~ fs/2
p = fftshift(p);
f = fftshift(f);
f(f>=fs/2) = f(f>=fs/2)-fs;

% 转换为dB（pwelch输出为PSD，因此是10倍log10）
pdb = 10*log10(p);

plot(f, pdb);
grid on;

在Cadence的计算器中，可以使用dft函数完成类似的操作，并且它也提供了分段、平均的选项来平滑输出结果，在此不再赘述。

CATALOG
    FEATURED TAGS
    MATLAB PLL Programing Math Others Flask Training SerDes ChatGPT
    FRIENDS