我们这届数字信号处理变成学院选修辣,相对的考试难度会下降一些。这篇blog记录一下我突击DSP时的一些重要考点。
由于这门课变成学院选修了,所以林老师其实删去了很多的考点,整体上变成了离散+阉割版信号与系统,再加上FFT的少许内容。
采样
这门课最重要的是理解“采样”,从连续的模拟量到离散的时间信号。这就引入了采样定理。证明我们就从略了,主要是理解。
“时域离散,频域延拓”。采样信号的频谱将会是原模拟信号的频谱沿频率轴,每隔一个采样角频率$\varOmega_s$重复一次,所以是周期性延拓。所以为了不出现频谱混叠,采样率$\varOmega_s$与模拟信号的截止频率$\varOmega_c$需要满足奈奎斯特定理,即:
因为采样后,频谱就被延拓了。(从另一个角度看,离散后的点实际上允许了更高频的信号去逼近它)。如果要从采样的离散点中恢复模拟信号,我们需要做插值(内插)。最简单的办法是离散点之间线性插值,但这个肯定是不对的。实际上,我们使用一个理想的低通滤波器,滤掉高频分量,即可。由时域卷积定理,这相当于在时域和一个$\mathrm{Sinc}$函数做卷积。只要满足奈奎斯特采样定理,这种恢复是无失真的。
很抱歉打断你的浏览,但我们要先明确一下,频率。
一个序列的周期是$T$,那么它的频率$f=1/T$,单位是Hz。角频率是$\varOmega=2\pi/T$,单位是rad/s。如果进行采样,我们会得到所谓的“数字角频率”$\omega$,举一个例子来说:
这里的$T$是采样周期,并不是信号的周期。它是采样率$f_s$的倒数,所以这个数字角频率$\omega$,可以表示为:
所以$\omega$在这里实际上是一个无量纲数,单位是rad。这点要注意区分。更进一步:
由采样定理,这个$\omega$是在$0\sim\pi$的。这个性质后面会用到。
系统的性质
在描述一个系统时,我们关注其是否满足周期性,线性,时不变性,稳定性,因果性。
周期性一般是对于离散后的正弦序列,只需计算其周期,如果为整数和有理分数,就是周期的。周期就是那个整数或有理分数的分子。
线性要满足:
所以带常数的系统就不是线性的了。
时不变性是:
所以当系统有伸缩,反转,以及内部的非线性时,就不满足时不变性。
稳定性,从最原本的定义上看,即有界的输入能保证有界的输出。特别地,对于离散时间的线性时不变系统,具有稳定性的充要条件是单位脉冲响应$h(n)$构成的级数绝对收敛。
因果性就很简单了,即系统的单位脉冲响应在$n<0$时都为0。这一点可能有些突兀,实际上,考虑:
用定义打开,就能看出来了。本质的意义是系统$n$时刻的输出序列,只取决于$n$时刻以及之前的输入序列。
离散时间傅里叶变换
DTFT只离散了时间,所以频谱实际上还是连续的。这点也能从上面的采样定理中看出,一般序列$x(n)$的DTFT定义为:
简记为:
由于在频谱上是连续的,逆变换时就是积分回去了,积分限是$-\pi\sim\pi$,原因见上面关于数字角频率的阐述。
考试只要掌握定义就好了,出了具体的序列,直接用等比数列梭哈:
计算序列$x(n)=R_N(n)$,求$x(n)$的DTFT:
这里$R_N(n)$是矩形序列,但不同于门函数,它直接定义的是从0开始到第$N$个序列点为1。这个表达比较方便,不然要讨论奇偶性。
阶跃序列并不是绝对可和的,所以和学信号与系统一样,求阶跃序列$u(n)$的DTFT需要特殊处理。首先我们要解决常序列的DTFT,即$x(n)=1$。它并不是绝对可和的,但在连续情形下:
注意这里的$\delta(n)$是0时为无穷大的冲激函数,并不是离散时$\delta(0)=1$的那个。我们对连续的频谱进行采样,频谱会周期延拓。这里我们取$T=1$,那么$\varOmega_s=2\pi/T=2\pi$,所以。对于离散后的$x(n)=1$这个不绝对可和的序列,就有了:
所以对于单位阶跃,我们有:
同时注意到:
所以我们得到:
给阶跃序列乘上一个$a^n,0<a<1$,得:
在学信号与系统时,我们多数都考虑输入信号为一个实函数。那时关于连续情形的傅里叶变换,只需记住一个“奇虚偶实”的结论,就可以了。即实函数的傅里叶变换,其频谱函数虚部是奇函数,实部为偶函数。
但在DSP这里需要进一步扩展这个点,书上写的,太乱了,叫人半懂不懂的:
红色部分就是之前奇虚偶实的来源,即考虑一个实数序列。那么它的DTFT只有共轭对称分量,共轭对称分量的虚部是奇函数,偶数是实函数。故奇虚偶实。
$x_e(n),x_o(n)$由于构造出来,天然有共轭对称和共轭反对称的性质。所以称$x_e(n)$也为偶序列(取对称,even之意),$x_o(n)$为奇序列(取反对称,odd之意)。
大多数时候,我们会关心离散时间的线性时不变系统的频率响应。输入序列$x(n)$通过单位脉冲响应$h(n)$,我们可以得到系统的零状态响应:
由时域卷积定理,可以导出:
例题我们合并到后面再给出~由于今年的考核不要求逆z变换(因为收敛域的问题),所以大概率会考单位脉冲响应的求解,那时候就需要逆DTFT了。
序列的Z变换
我们作共形映射$z=e^{j\omega}$,DTFT的基被映成了一个复平面上的单位圆。扩展$z$的范围,就得到了一个洛朗级数:
在考核里我们重点关注用$z$变换求解差分方程,由于一般输入$x(n)$都是因果的,所以单边还是双边$z$变换结果是相同的。
在过去,我们曾引入连续的拉普拉斯变换,然后再到$z$变换,再到收敛域,因果性,稳定性。今天我们可以换一个角度。
考虑因果系统的充要条件,是单位脉冲响应$h(n)$在$n<0$时等于0,所以考虑:
那么$z$的收敛于一定包含无穷远处。所以收敛于是某个圆的外部。
而稳定系统的充要条件是$h(n)$绝对可和,这同时允许$h(n)$的DTFT存在。所以,当系统稳定时,收敛域一定包含单位圆。
综合以上两点,对于一个因果稳定系统,其所有极点一定在单位圆内。我们也可以通过上面的结论,从极点分布和收敛域上,分析一个系统的因果性和稳定性。
引入零点和极点表示还有一个好处,它可以显式的给出获得系统频率响应的几何方法。结论是:系统函数的极点位置主要影响系统幅频响应的峰值和尖锐程度,零点位置主要影响系统幅频响应的谷值和谷深。这个结论其实从一定程度上看是显然的。
最后,我们强调收敛域对于$z$变换的意义,同时引出一道可能考察的类型。
在$z$变换中,我们不要求序列绝对可和,所以,就可以真的直接级数求和的梭哈了,吗?考虑$x(n)=a^n u(n)$的$z$变换及其收敛域:
再考虑$x(n)=-a^n u(-n-1)$的,会发现:
此时的收敛域,是$\left| z \right|<\left| a \right|$,所以,不同的序列,$z$变换的式子是一样的。如果不加约束,直接给出一个$X(z)$,那只要不包含极点的都可以是它的收敛域,所以需要分情况讨论,指定收敛域。
例如一个线性时不变系统由下面的差分方程描述:
确定其系统函数$H(z)$,给出其收敛域,作出其零极点分布图,再分析该系统的稳定性。
我们直接计算:
分子显示零点为$z=0,z=0.4$,极点为$z=0.5,z=2$。
这里零极点图就不画了。(零点是圈,极点是叉)
所以收敛域有三种情况:
当$\mathcal{R}_1$时,系统由因果序列和反因果序列共同构成,此时收敛域包含单位圆,系统是稳定的;当$\mathcal{R}_2$时,系统只剩因果序列,因为收敛域包含无穷远,由于没有单位圆,系统不稳定;当$\mathcal{R}_3$时,系统是反因果序列,不包含单位圆,不稳定。
离散傅里叶变换
闪回到离散时间傅里叶变换,我们知道:
频谱函数还仍然是连续的,如果进一步,我们对频谱函数进行$N$点等间隔采样,我们也就实现了频域上的离散。到这一步,我们所构想的频域才能真正在计算机上实现。
我们考虑这个离散的过程,由于$\omega$是属于0到$2\pi$的,所以一个自然的分度是$2\pi/N$,(这个$N$也恰好是序列的采样点数,我们后面会讨论为什么这两个数值一般的时候是相等的。)而原有序列$n$的索引仍然可以以指数的形式存在。同时第$k$个频率可以表示为$\omega_k=\frac{2\pi}{N}k$。
我们用新的记号来表示等间隔采样后的$e^{-j\omega n}$:$W_{N}^{kn}=e^{-j\frac{2\pi}{N}kn}$。沿着$k$建立索引,就得到了离散傅里叶变换:
这里求和只作用在了主值区间中。形象表示如下图所示:
实际上,$W_{N}^{kn}$可以写成矩阵形式:
这个矩阵的行列式是一个范德蒙德行列式,由于$x_i,x_j$互不相同,所以其行列式不等于0。所以它是可逆且唯一地确定了$X(k)$与$x(n)$之间的映射。
DFT也存在共轭对称性,关于这一性质,我们在DTFT中有详尽的讨论。在DFT中,如上面体现DFT的插图所示,我们的对称性从关于原点转移到了关于$N/2$,其余性质都得以保留。我们以实序列的共轭对称性为例:
如果$x(n)$是长度为$N$的实序列,那么像在DTFT中讨论的那样,$X(k)$只有共轭对称分量,所谓共轭对称反映在:
所以我们只需求得前一半数目的$X(k)$,就能得到另外一半的了。特别地,如果$x(n)$在实序列的基础上还是个偶对称的,那么重复利用在DTFT中发现的性质,其共轭反对称分量$x_o(n)$为零,则其DFT不存在虚数部分。所以$X^*\left( N-k \right)=X^\left( N-k \right)$,于是$X(k)$也就有了实偶对称:
同理,如果序列是实奇对称。那么不存在$x_e(n)$,其DFT只有虚数部分。所以可以将共轭直接换成负号,此时$X(k)$呈纯虚奇对称:
更重要的是,我们可以用DFT进行一些更一般的频谱分析。而不是拘泥于先前,只能计算一些特定序列的级数和。我们现在着重考虑用DFT分析连续时间信号,我们先着重分析这张图的横坐标:
最上面是模拟信号的频谱$X_a(j\varOmega)$,此时的横坐标是角频率$\varOmega=2\pi/T$,其截止频率是$\varOmega_c$。中间是进行采样,得到的$X(e^{j\omega})$,我们先前讨论过,它其实是角频率被采样频率$f_s$归一化的结果$\omega =\frac{\varOmega}{f_s}$,由于采样频率$f_s$一般大于等于2倍的$f_c$,所以我们看到$X(e^{j\omega})$主值区间截止于$T\varOmega_c<\pi$。这里的$T$是采样周期,即$1/f_s$。
接下来,$X(e^{j\omega})$被直接进行了抽取,间隔为$2\pi/N$,我们现在关注在$X(k)$中第$k$个频率对应在$X_a(j\varOmega)$的频率是多少。这需要分情况讨论,因为$X_a(j\varOmega)$没有经过周期延拓。
当$0\le k\le N/2-1$时,我们可以得到:
当$N/2\le k\le N-1$时,我们需要减去一个周期:
更进一步,我们可以把大于$N/2$的那些$k$,平移到负半轴上。这样重排一下就变得对称了,同时,$k$与$\varOmega_k$也有了一个不间断的映射:
明确了这一关系,我们再来看DFT在近似频谱分析时的问题,以及解决方法。
①混叠现象:
我们知道,时域有限的信号,频谱无限宽;频谱有限宽的信号,持续时间无限长。所以,持续时间有限且带限的信号,其实不存在。如果信号的频谱很宽,直接进行采样再DFT,需要很高的采样率才能不导致混叠。所以另一个办法是采样前利用一个模拟的低通滤波器进行预处理,限制其上限频率,以减少混叠程度。
②截断效应:
如果一个信号持续时间很长,我们只能对其进行截断,得到一个有限长的序列。这个过程形象的叫作时域加窗,即我们对信号乘上了一个窗函数:
时域相乘相当于频域卷积,如果我们不加任何处理,我们其实选取的就是一个矩形窗。它的傅里叶变换是一个$\mathrm{Sinc}$函数,其幅度值会呈现“主瓣”和“副瓣”,就像:
我们以一个余弦的频谱为例,体会一下其影响:
我们得以观察到两个现象,一是频谱泄露,由于一个$\mathrm{Sinc}$卷积的作用,冲激函数附近产生了较大的扩展,这种现象有时会淹没掉频率相近但微弱的信号。二,谱间干扰,这种扩展本身也会造成不可避免的混叠,有时他们会混叠出一个峰值,可能会被误认为是另一个信号的谱线,造成误判。
这种截断效应,可以通过增加截断长度来缓解。还可以通过使用主瓣更宽,副瓣更窄的窗函数,这些在后面会介绍。
③栅栏效应:
由于DFT是DTFT的等间距采样,我们相当于透过一个栅栏在观察频谱。所以有可能一些峰值和谷点就被挡住了。这称为栅栏效应,减少它的方法就是对原序列补零,再作DFT。等价于在频率上增加采样的点数。
跟这个有些类似的概念是:提高信号的频谱分辨率。看似,我们增加采样点数,可以提高频率的分辨率,这是一种误解。因为频率的分辨率,在模拟信号经过采样,截断后,已经确定了。为$f_s/N$,补零只能提高频谱显示的分辨率,并不能真的恢复那些已经被损失的频率信息。
举个例子就是,你是无法通过补零,来把两个频率相近的小信号和大信号的频谱分开的。提高信号频率分辨率的方法只有从根本上提高采样频率或增加截取长度。
现在,我们结合上面的现象,给出一个利用DFT进行信号分析时的参数选择。在大多数情况下,信号的截止频率$f_c$是已知的(如果它带限,我们就可知;不带限用抗混叠滤波器过一遍,也已知了),同时我们会对信号的频率分辨率$\varDelta f$提出要求。我们只考虑矩形窗加窗的情况,由上图我们可知主瓣的有效宽度是$2\pi/N$。
所以,我们要求主瓣的宽度要小于频率分辨率,这样才不至于出现谱间干扰,于是:
采样率,很自然的要满足$f_s\ge2f_c$;信号的持续时间$T_p$(这是相对采样前说的)应该满足分辨率的要求,和上面推理$N$的式子是类似的:
这个考试可能会出分析和计算题。
线性卷积和循环卷积
我们很久之前就听过:“可以用傅里叶变换加速卷积运算”,那只是一个理论上的剪影,我们并没有从工程上思考过这件事。
首先我们要明确,我们熟悉的卷积运算:
是线性卷积,使用DFT来实现序列的线性卷积,我们还需要作一些准备工作。
首先,我们定义一种移位运算,它很直接,看图一下就能看懂:
这种运算指的是:
在这个基础上,扩展了一种循环卷积,这种循环卷积可以满足在离散序列下“时域卷积等价于频域相乘”。
证明这里就不管了,我们重点关注使用DFT来计算线性卷积。可以证明,对于长为$N$的序列$x(n)$和长为$M$的序列$y(n)$,其循环卷积点数$L$需要满足$L\ge N+M-1$,此时循环卷积才等于线性卷积。
我们用例题来记忆上述过程,已知:
计算两序列线性卷积$y_L(n)$,7点循环卷积$y_C(n)$,并作出由DFT计算线性卷积的框图,标出每一步所需的点数。
线性卷积是好计算的,我们直接用不进位乘法就可以得到:
此时由于满足$7\ge 4+4-1$,所以循环卷积的结果和线性卷积是一样的。框图即为:
其中$M=N=3, L=7$。
下面我们给出一种更一般的方法来计算循环卷积,它是说先算出对位相乘相加序列,然后根据$L$修正。若$L$为4,则把4之后的序列搬移到前面,两序列求和即得最终结果。
例如$x_1(n)=\{-3,2,4,1\},x_2(n)=\{1,1,1\}$,作$x_1\left( n \right) \circledast x_2\left( n \right) $,$L=4$。
我们得到:
快速傅里叶变换
考试中要求作出FFT的蝶形计算流图,今年只用记住时域抽取基2FFT算法即可。
这里,要记住三个事情:
输入$x(n)$的顺序:使用“逆序二进制加法”,例如$N=8$时,从000开始,先最高位加1,得100(4),然后继续最高位加1,向次高位进位,得010(2),然后是110,6。所以顺序即为$x(0),x(4),x(2),x(6)…$
第$L$级旋转因子$W$的填写:按照$W_{N}^{J\cdot 2^{M-L}}$,其中$J=0,1,2…(2^{L-1}-1),M=log_2N$。例如第一级,$J=0$,所以都是$W_8^0$,到了第二级,$J=0,1$,所以$W_8^0,W_8^2$,以此类推……
与DFT相比的计算量:
如果使用传统的矩阵相乘,需要有$N(N-1)$次复数加法,$N^2$次复数乘法;用FFT计算,如蝶形图所示,乘法次数为$\frac{N}{2}\log _2N$,加法次数为$N\log _2N$。
实际上,这个蝶形单元并不是一个正式的信号流图。所以,照做吧,死记硬背。
End
这才仅仅是一半,另一边是和滤波器设计相关的,记在下一篇blog了。