生物医学信号处理之数字信号处理基础

数字信号处理基础-离散傅里叶变换（DFT）

本文主要参考书籍：
[1] 饶妮妮，李凌. 生物医学信号处理[M]，成都，电子科技大学出版社，2005.7.
[2] 姚天任，孙洪. 现代数字信号处理 （第二版）[M]. 武汉华中科技大学出版社，2018.8.

前言

大多数理工科高等学校分别为本科生和研究生设置课程“数字信号处理”和“现代数字信号处理”。粗略地说，其课程内容是按照所要处理的信号和处理方法来划分的。具体说，前者针对离散时间确定性信号，而后者针对离散时间随机信号。因此，前者以数字信号处理学科的两大理论支柱，即离散傅里叶变换和数字滤波器作为核心内容；而后者以自适应滤波器（维纳滤波器和卡尔曼滤波器可以看成是其基础）、功率谱估计、小波分析（一种最典型、最有力的时频分析方法）、同态信号处理、高阶谱分析，以及神经网络信号处理为主要内容。

    这篇博文主要是介绍前者-数字信号处理，后者-现代数字信号处理将在后续文章中整理分享。

    对采集到的信号进行处理，除了传统的时域分析之外，各种各样的变换发挥了重要作用，包括从最熟悉的傅里叶（Fourier）变换到现在的小波（Wavelet）变换以及主成分分析（Principal Component Analysis）、独立成分分析（Independent Component Analysis）和稀疏成分分析（Sparse Component Analysis）。每一种变换都有其独特的视野，为信号的分析处理提供了不同的思路。这里我们将介绍最基本的一种变换，即在线性时不变系统（LTI，Linear Time Invariance）中广泛使用的傅里叶变换以及频谱分析，本章还是以离散数据为主，介绍离散傅里叶变换（DFT）的有关知识。

1 傅里叶变换

    傅里叶分析方法的处理对象包括连续时间信号和离散时间信号，在近代欧拉、伯努利、傅里叶、狄里赫利等学者努力完善了对连续时间信号的傅里叶分析方法。用于处理离散数据以产生数值近似的有关内插、积分和微分等方面的公式早在17 世纪的牛顿时代就被研究过，从事时间序列的研究曾吸引了18、19 世纪包括高斯在内的许多著名科学家，从而为离散傅里叶变换提供了数学基
础。
    在20 世纪60 年代中期，库利（Cooley）和图基（Tukey）各自独立发表了一篇论文，也就是快速傅里叶变换算法（FFT）。FFT 是非常高效的算法，使得计算变换所需要的时间减少了几个数量级。由于计算机速度的迅速提高，越来越多的连续时间信号被离散化，然后用计算机进行处理。

    利用“三角函数和”的概念来描述周期性过程，三
角函数和也即是成谐波关系的正弦和余弦或周期复指数函数的和。如果一个LTI 系统的输入可以表示为周期复指数的线性组合，则输出也一定能表示成这种形式，并且输出线性组合中的加权系数与输入中对应的系数有关。
    如图2-1 所示，x(n)表示输入或者激励， y(n)表示系统输出或者响应，$H(e^{jw})$表示系统单位脉冲响应h(n)的频率响应。

    在研究LTI 系统时，复指数信号的重要性就体现在图2-1 中：一个LTI 系统对复指数信号的响应也同样是一个复指数信号，不同的只是乘了一个复振幅因子$H(e^{j{w}_k})$，频率并没有发生变化，由于是复数因子，就有了幅度和相位或者实部和虚部的变化。

连续和离散时间信号的傅里叶级数和傅里叶变换表达式如下：

    从上表可发现傅里叶变换在LTI 系统分析中的思想，就是把一个无论多复杂的输入信号分解成复指数信号的线性组合，那么系统的输出也能通过如图2-1 所示的关系表达成相同复指数信号的线性组合，并且在输出中的每一个频率的复指数函数上乘以系统在那个频率的频率响应值。

$\tilde{x}(t)$表示时域连续周期信号
$x(t)$表示时域连续非周期信号
$\tilde{x}(n)$表示时域离散周期信号
$x(n)$表示时域离散非周期信号

系数 { a k } { \{a_k\}} {ak​} 称为连续周期信号$\tilde{x}(t)$的傅里叶级数系数或频谱系数或线谱等,是离散非周期的；
系数 { a ~ k } {\{\widetilde{a}_k\}} {a k​} 称为离散周期信号$\tilde{x}(n)$的傅里叶级数系数或频谱系数或线谱等,是离散周期的；
$X(jw)$称为连续非周期信号x(t)的频谱，是连续非周期的；
$X(e^{jw})$也称为离散非周期信号x(n)的频谱，是连续周期的。

    作为线性组合所取的形式从求和过渡到积分，就是利用傅里叶的思想，一个非周期信号可以看成是周期无限长的周期信号，当周期增加时，基频$w_0$ 越小，成谐波关系的各分量在频率上越来越近，当周期变得无穷大时，离散的线谱就形成了一个连续谱，也就从求和变成了积分。

    从表2-1 中时域和频域的关系还能得到如下规律：时域的离散必然导致频域的周期化；频域的离散必然导致时域的周期化。简单地说，就是一个域离散必然另外一个域是周期的；相反，如果一个域连续必然另外一个域是非周期的。掌握了这个规律，我们很快就能判断出一个信号在频域的表现形式。
   
傅里叶级数对应的频域是离散的，时域是周期的；傅里叶变换对应的频域连续的，时域是非周期的。

    由于计算机无法处理连续信号，因此需要一种在时域和频域上都离散、非周期的一对傅里叶变换对，这就是离散傅里叶变换，简称DFT。离散傅里叶变换的到处方法有多种，比较方便同时物理意义清晰的是从离散时间傅里叶变换（DTFT）和从离散傅里叶级数（DFS）入手。
离散时间变换DTFT为：
$$X(e^{jw})=\sum _{n=-\infty }^{+\infty }x(n)e^{-jwn}\text{\hspace{1em}(1)}$$
时域是离散的，非周期的，但频域是连续、周期的，对连续变量$w$均匀采样，也就是对单位圆进行$N$等分，取一个周期的结果即得：
$$X(e^{jw}){|}_{w=-\frac{2\pi }{N}k}=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk},(0\le k\le N-1)\text{\hspace{1em}(2)}$$这样频谱变量由连续量$w$变成了离散变量。

从DFS到DFT更加明显，DFS对应的时域和频域都是离散周期信号，可以在这两个与分别取他们的主值，也就是限定在一个周期内，这样就得到了DFT的变换对。具体给出如下：
$$X(k)=DFT[x(n)]=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk},(0\le k\le N-1)\text{\hspace{1em}(3)}$$

$$x(n)=IDFT[X(k)]=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk},(0\le n\le N-1)\text{\hspace{1em}(4)}$$

称（3）式为正变换，（4）称为反变换。

注意：这一对变换对中信号x(n)的长度为N，它的频谱X(k)点长也为N，则x(n)和X(k)具有唯一的映射对应关系。也有可能给一个7 点的时间序列，求该信号的4 点DFT 或者10 点DFT，明显前者求得的频谱不能真实反映出信号，而后者不存在混叠现象，能反映出信号的频谱。因此在求DFT 时要注意到底是求多少点长的DFT，只有时域样点数小于或等于频域样点数，频谱才是真实的反映。
下面给个实用易懂的例题：

2 傅里叶变换的性质

    离散傅里叶变换是有限长序列的Z 变换在单位圆上抽样的结果，因而很多性质和序列
    傅里叶变换性质类似，但由于DFT 又可以从离散傅里叶级数引入，因而隐含了周期性，所
以又有些不同。

帕塞瓦尔（Parseval）定理

3 频域分析和图谱表示

    离散傅里叶变换作为傅里叶变换的一种近似而得到广泛应用，它的快速算法保证了DFT 在实时信号处理中的应用。下面介绍频谱分析中常用的几种方法。

3.1 幅度谱

设N 点长序列$x(n)$的DFT结果是$X(k)$，它是离散的复序列，可以用下式表示：
$$X(k)=|X(k)|e^{j\angle X(k)}\text{\hspace{1em}(3-1)}$$
代入傅里叶变换的反变换公式有：
$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}|X(k)|e^{j[\frac{2\pi }{N}nk+\angle X(k)]}\text{\hspace{1em}(3-2)}$$
离散傅里叶变换的模$|X(k)|=\sqrt{[Re(X(k)){]}^2+[Im(X(k)){]}^2}$ ，表示信号$x(n)$的各复指数信号的频率分量$（w_k=\frac{2\pi }{N}k,0\le k\le N-1）$的相对大小。
例如，如果信号$x(n)$在k = 0 附近小范围以外$X(k)=0$，那么$x(n)$所呈现的仅是相当低的频率。
如果幅度谱是偶对称的，往往只需要画一半即可。
画幅度谱时，采用对数坐标也是很常用的，即幅度大小用$20lo{g}_{10}|X(k)|$来代替，这时纵坐标的单位就是分贝（dB），0 分贝对应模等于1，20 分贝就对应10 倍的增益，−20 分
贝对应于衰减0.1，等等。
例题附上，结合矩形窗DFT的例题有助于理解：

3.2 相位谱

$\angle X(k)$表示相位角， $$\angle X(k)=t{g}^{-1}(\frac{Im[X(k)]}{Re[X(k)]})\text{\hspace{1em}(3-3)}$$
它的大小不会影响各复指数频率分量的大小，但能提供这些频率的初始相位信息。$\angle X(k)$对信号$x(n)$的性质有着显著的影响，因此一般包含了信号的大量信息，用相同的幅度谱和不同的相位谱得到的信号完全不同。
如果序列$x(n)$是实序列，$\angle X(k)=-\angle X(N-k)$，即相位谱是奇对称的。
【例2-7】画出例2-6 的相位谱并解释该图。
解：因为：

3.3 功率谱

    信号$x(n)$的离散傅里叶变换$X(k)$一般是一个复数，$X(k)$与其共轭$X^{\ast }(k)$之积称为自功率谱，简称自谱或功率谱，或功率谱密度函数，其表示为：$$P(k)=\frac{1}{N}X(k)X^{\ast }(k)=\frac{1}{N}|X(k){|}^2,0\le k\le N-1\text{\hspace{1em}(3-4)}$$

    系数$1/N$是为了满足能量定理而进行的调整。P(k)反映的是信号的功率密度，在图形上与幅度谱类似，只是大小不同，功率谱不含相位信息，所以由功率谱不能恢复原始信号，为存在多义性。

4 频域分辨率

   利用离散傅里叶变换进行频谱分析时会引起误差，常见的就是由于离散傅里叶变换本身采样和截断过程所引起的混叠、泄漏、栅栏现象等。

   如果要分析的信号是连续信号，就必须先采样，当采样频率比信号最高频率的两倍要小，即fs＜2Fh 时，就会发生混叠现象，可以提高采样率来避免混叠现象。

   如果要分析的信号是周期连续信号，就必须对该信号截取一段来进行分析，即加了一个窗，便会发生泄漏现象，这时离散傅里叶变换的结果似乎是将原信号频谱扩展到整个频率范围，使得无法反映真正的频谱，要想减少泄漏可以通过加不同的窗函数来截取信号。

   离散傅里叶变换是对离散时间傅里叶变换的采样，它只给出频谱在离散点$（w_k=\frac{2\pi }{N}k)$上的值，而无法反映这些点之间的频谱内容，这就是栅栏现象，改善栅栏效应的一种方法是信号后面补若干个零，通过补零来调整坐标上$w_k=\frac{2\pi }{N}k$ 的位置。

   在时域中，数字信号的时间分辨率由采样间隔Δt 决定。在频域中，频域分辨率$\Delta f$（两相邻谱线间的频率差值）由采样频率fs 和采样点数N 决定：
$$\Delta f=\frac{f_s}{N}\text{\hspace{1em}(3-5)}$$

式中，$f_s$ 由采样定理决定，因此，要提高频域分辨率，就得增加采样点数N。如果数据量略有不足，传统方法是在数据尾部填0 来解决，称为高密度频谱（The High Density Spectrum）。但是补零并不能提高频域分辨率，我们认为填入适当的现有数据会更好，称为高分辨率频谱（The High Resolution Spectrum）。
   一般来说，功率谱、幅度谱、位相谱等皆为离散值而非连续值，把上述谱绘成连续谱线图是不恰当的，是误解，在由频率分辨率决定的两条谱线间的值是不能确定的或不存在的，如果将谱线绘成连续曲线，则会误解为两条谱线间的值是确定的。假如频率分辨率为$f_0=1Hz$，而原始信号中又确有如1.3Hz 或1.7Hz 的频率成分存在，在谱线中是表现不出来的。
例题贴上：


通过上述讨论，利用FFT 做谱分析时各参数的选择为：
1．采样频率应满足采样定理：
$$采样频率f_s\ge s{F}_h\text{\hspace{1em}(3-6)}$$
或者 $$采样间隔T\le \frac{1}{2}{F}_h\text{\hspace{1em}(3-7)}$$
往往一段信号的频谱是在整个频率范围内，因而要先通过一个带通滤波器，来限制信号的最高频率。

2. 离散傅里叶变换的点长N 为
$$N=\frac{f_s}{\Delta f}\ge \frac{2F_h}{\Delta f}\text{\hspace{1em}(3-8)}$$
一般，为了FFT 计算快速，N 都尽量取成2 的幂次。

3．采集信号的持续时间为
$$p_t=NT\text{\hspace{1em}(3-9)}$$

5 数字里滤波器

理想的数字滤波器幅度谱如图2-8 所示，这些频率特性都是以2π为周期的连续函数，当单位脉冲响应h(n)是实数序列时，幅度谱周期偶对称，相位谱周期奇对称，因而只需要给出一半的频谱图即可。

数字频域滤波可以用硬件和软件的方法实现，滤波器的设计方法多种多样，大致分为IIR 滤波器的设计和FIR 滤波器的设计，前者主要利用传统的模拟滤波器设计方法，后者多采用窗函数和频率取样设计法.
下面结合Matlab 举例说明数字滤波器的设计实现。
【例2-11】选择一种合适的窗设计一个FIR 数字低通滤波器，满足下列要求：
$$w_p=0.2\pi ,A_p=0.25dB,w_s=0.3\pi ,A_s=50dB$$
并画出滤波器的单位脉冲响应和该滤波器的幅度响应。
解：下面用Matlab 语句实现该滤波器的设计以及画出相应的图形。

wp = 0.2*pi; ws = 0.3*pi;  %给出通带频率和阻带频率
tr_width = ws–wp; %求过渡带宽度
% 50dB s QA = ，hamming window即可满足该条件，查表求得窗长度
M = ceil(6.6*pi/tr_width) ;
n = [0:1:M−1];
wc = (ws+wp)/2; %求截止频率
b= fir1(M,wc/pi); %求FIR 低通滤波器的系数，默认就是hamming window
h = b(1:end−1);
[hh,w] = freqz(h,[1],'whole'); %求滤波器的频率响应
hhh = hh(1:255);ww = w(1:255); %由于对称性，画一半图即可
%画图
subplot(1,2,1); stem(n,h);title('实际脉冲响应')
axis([0 M−1 −0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(1,2,2); plot(ww/pi,20*log10(abs(hhh)));title('幅度响应(单位：dB)');grid
axis([0 1 −100 10]); xlabel('频率(单位：π)'); ylabel('分贝')
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])
set(gca,'YTickMode','manual','YTick',[−50,0])
%图形如2-10 所示。

【例2-12】最常碰到的信号处理任务是平滑数据以抑制高频噪声。求几个数据点的平均值是减弱高频噪声的一种简单方法，这种滤波器称为平滑滤波器或中值滤波器。下面给出一段1000 点的脑电数据，用不同长度的数据点进行平滑，观察其滤波后的效果。
解：在Matlab 中有直接的函数进行平滑：$Y=MEDFILT1(X，N)$，如果没有给出N 的值，则默认N= 3；当N 是奇数时，Y 是$X(k-(N-1)/2:k+(N-1)/2)$的平均；当N 是偶数时，Y 是$X(k-N/2:k+N/2-1)$的平均。这里给不同的N 看其滤波后的效果：

n = 0:999;
y1 = medfilt1(x, 15);
y2 = medfilt1(x, 30);
y3 = medfilt1(x, 45);
subplot(2, 2, 1); plot(n,x); legend ('original signal');
subplot(2, 2, 2); plot(n,y1); title('N=15'); legend ('filtered signal');
subplot(2, 2, 3); plot(n,y2); title('N=30'); legend ('filtered signal');
subplot(2, 2, 4); plot(n,y3); title('N=45'); legend ('filtered signal');
%如图2-11 所示，N 越大，曲线越来越光滑。

Matlab 中有多个滤波器设计和实现的函数，例如：
FIR 滤波器的设计：fir1、fir2、firls、remez、fircls、fircls1、firrcos、sgolay 等。
IIR 滤波器的设计：butter、cheby1、cheby2、ellip、yulewalk 等。
滤波实现：filter、medfilt1、filtfilt、freqz、freqs、latcfilt 等。
具体使用方法可以用Matlab 的帮助系统查询。

作业完成，下一次梳理随机信号基础，期待~

本文主要参考书籍：
[1] 饶妮妮，李凌. 生物医学信号处理[M]，成都，电子科技大学出版社，2005.7.
[2] 姚天任，孙洪. 现代数字信号处理 （第二版）[M]. 武汉华中科技大学出版社，2018.8.