通信系统的计算机模拟第六讲

通信系统的计算机模拟第六讲1第一页，共五十七页，2022年，8月28日第五讲回顾带通信号低通仿真模型复包络：频域，同相和正交分量的谱能量与功率随机带通信号的正交模型随机复包络信号的若干结论信噪比与系统的低通仿真模型线性时不变系统多载波信号非线性与时变系统2第二页，共五十七页，2022年，8月28日第五章滤波器模型与仿真方法滤波器是组成通信系统的许多子系统中的重要部分这些滤波器中的许多是模拟的，为了便于仿真，必须将它们映射为合适的等价数字滤波器近似并且都会在仿真结果中引入误差根据定义，滤波器具有频率选择性，并且具有有限时间冲激响应或者无限时间冲激响应。由于滤波器具有频率选择性，它们引入了记忆???，而由于这种记忆，需要滤波器过去的输入和/或输出来计算给定时刻的滤波器输出。因此，滤波器需要存储，而存储和恢复采样值会大大地增加仿真程序的计算量，反过来，这又增加了仿真程序的运行时间，因此，我们寻求较好的仿真结构以获得减少计算量的算法。不是滤波器设计3第三页，共五十七页，2022年，8月28日5.1概述数字滤波器的经典综合方法通常都是基于模拟滤波器原型之上的仿真模型所要求的数字滤波器的设计，从模拟滤波器的S（拉普拉斯）域传递函数着手，来寻求等价的数字滤波器，等价于模拟原型的数字滤波器。基本方法是基于时域判据或频域判据来衡量其等价性。通过要求数字滤波器的输出与作为设计基础的模拟原型的采样输出相匹配，来运用时域判据进行滤波器的设计。“冲激不变法”是使数字滤波器的冲激响应与模拟原型的冲激响应的采样相匹配的一种设计方法对于阶跃不变数字滤波器，数字滤波器的阶跃响应则与对应的模拟原型滤波器阶跃响应的采样相匹配。模拟原型和得到的数字滤波器在时域上是等价的，它们在频域上也近似等价。（采样频率混叠）4第四页，共五十七页，2022年，8月28日双线性z变换最常用方法：双线性z变换：模拟原型（双线性z变换）数字等价模型它允许对频率响应上的点进行匹配，使得模拟原型和所得数字滤波器在某些指定频率值上具有相同的频率响应（包括幅度和相位）。双线性变换法以引入非线性频率畸变（Warping）为代价消除了混叠误差。5第五页，共五十七页，2022年，8月28日IIRandFIR基于模拟原型的综合方法给出无限时间冲激响应（infinite-durationimpulseresponse,IIR）数字滤波器。原因：模拟滤波器的冲激响应会渐近趋于零，严格地讲，这种冲激响应在持续时间上是无限长的。为获得用于仿真程序的有限时间冲激响应（finite-durationimpulseresponse,FIR）数字滤波器，由标准IIR设计方法所产生的冲激响应可能需要截短（truncate）。通过在FIR滤波器的冲激响应中包含足够多的项数，由截短引起的误差可以降低到一个允许的水平下。6第六页，共五十七页，2022年，8月28日Pure数字滤波器数字滤波器的一个重要特性是，在没有对应的模拟滤波器时也可设计数字滤波器。允许近似给定的幅度响应，但要保持一个完全线性的相位响应，这种滤波器为FIR滤波器，可通过横向延时线（TransversalDelayLine,TDL）结构实现。这种滤波器有多种设计方法，最基本的方法是将以采样频率为周期的期望幅度响应展开成为傅里叶级数，所得傅里叶级数的系数就确定了数字滤波器的冲激响应。这是频率采样的一个例子——因为在各频率点上对期望频率响应进行“采样”，对这些频率采样点进行傅里叶反变换（IFFT）之后，就获得了滤波器的冲激响应。将滤波器的输入信号与冲激响应进行卷积运算，实现了仿真滤波器模型并且获得了滤波器的输出信号。7第七页，共五十七页，2022年，8月28日数字滤波器分类8第八页，共五十七页，2022年，8月28日

5.2IIR与FIR滤波器数字滤波器通常根据冲激响应的持续时间分为IIR滤波器和FIR滤波器。和冲激响应分类紧密联系的是滤波器的实现或结构。9第九页，共五十七页，2022年，8月28日5.2.1IIR滤波器数字滤波器：计算过去的N个输出采样点y[n-k](1kN）、当前的输入采样点x[n]，以及过去的N个输入采样点x[n-k]](1kN)的加权和作为当前的输出采样点。任何权值ak或bk(k0)具有非零值，处理器就会有记忆，因而具有频率选择性，我们称之为“滤波器”。在时变系统中，其中的一个或多个权值还会是序数n的函数。对（5-1）式两边进行Z变换，可以得到传递函数H（z）。记住z变换是一种线性运算，和的z变换等于z变换的和，k个采样周期的延时等价于乘以因子z-k。（5-1）10第十页，共五十七页，2022年，8月28日由上式得传递函数为此即线性时不变滤波器的传递函数的一般形式（5-2）（5-3）11第十一页，共五十七页，2022年，8月28日滤波器的冲激响应以及频率响应冲激响应（记为h[n]）是传递函数H（z）的逆变换，用exp(j2πfT)替代传递函数中的变量，就可得频率响应。数字滤波器的冲激响应记作h[n]，只要令x（n)=δ(n)就可得到它。

（5-4）（5-5）12第十二页，共五十七页，2022年，8月28日WhyIIR由于式（5-1）的递归性(y[n]是y[n-1]的函数），冲激响应通常是无限长的，IIR滤波器。冲激响应h[n]是序数n的离散函数，而频率响应是连续变量f的连续函数。13第十三页，共五十七页，2022年，8月28日

5.2.2FIR滤波器对于k≥1有ak=0,则得到FIR滤波器。FIR数字滤波器的传递函数如下：其对应的冲激响应下如：仅在取值范围0≤n≤N内，冲激响应具有非零值，最多有个N+1非零项，它是有限长的。（5-6）（5-7）（5-8）14第十四页，共五十七页，2022年，8月28日5.2.3综合与仿真在仿真程序中使用滤波器涉及两种截然不同的运算操作综合：给定了滤波要求，并确定了满足这些要求的滤波器的传递函数H（z），这就建立了仿真模型。综合运算所得结果通常表示成两个向量：一个包含分母系数，另一个则包含分子系数。这两个向量定义了传递函数式，以及给定滤波器的输入产生滤波器输出的算法。即使采用复数算法，综合操作的计算量通常也不大，因为滤波器的综合运算只进行一次，因而发生在仿真程序主循环的外部。第二种运算操作涉及在每一个仿真时间步，计算滤波器的输出。这种运算在蒙特卡罗仿真程序中要重复进行数百万次，甚至数十亿次，如果想在合理的仿真运行时间内完成，要尽可能地减小这种运算的计算量。15第十五页，共五十七页，2022年，8月28日5.3IRR与FIR滤波器的实现Target:最小化计算量仿真时间短5.3.1直接II型与转置直接II型的实现IIR数字滤波器的一种有效方法是采用转置直接Ⅱ型结构。转置直接Ⅱ型滤波器结构的信号流图，及推导出这种结构的直接Ⅱ型结构的信号流图。16第十六页，共五十七页，2022年，8月28日17第十七页，共五十七页，2022年，8月28日产生转置直接Ⅱ型滤波器的规则1．重画原有（直接Ⅱ型）信号流图，保持结构（所有链路的相对位置不变）。2．将每次链路中的信号流反向。3．对新链路赋上与原有链路相同的运算（乘以一个常数、时延等）。4．如有需要，翻转（从左到右）新构造的信号流图，使得输入/输出信号流的方向与原始信号流图一致（注：信号流的方向通常是从左向右）。新获得的信号流图叫做转置直接Ⅱ型（DFⅡ）结构，其传递函数与原始信号流图相同。18第十八页，共五十七页，2022年，8月28日转置直接Ⅱ型的优势5-2所示的四阶转置直接Ⅱ型滤波器。（可以很方便地扩展到高阶滤波器）。给定输入采样，要计算输出信号第一步，计算状态变量。对于四阶滤波器，在我们的公式中有五个状态变量。这五个状态变量的计算式如下：（5-9）（5-10）（5-11）（5-12）（5-13）19第十九页，共五十七页，2022年，8月28日直接对应着Matlab程序20第二十页，共五十七页，2022年，8月28日x和y分别代表滤波器的当前输入和当前输出。式（5-9）到式（5-13）采样步长的延时，不需要存储和恢复这是MATLAB子程序filter要基于转置直接Ⅱ型结构的原因。在首次进入仿真循环之前必须进行初始化，输出一个瞬态响应。这段时间通常称为“稳定时间”，其大小为滤波器带宽倒数的若干倍。21第二十一页，共五十七页，2022年，8月28日矩阵形式数字滤波器的状态方程通常使用矩阵形式来表达W[n]和W[n-1]是的k×1的列向量，分别表示当前及过去的状态变量。Fc和Fd是k×k的系数矩阵，B是k×1的列向量，用来将输入x[n]耦合到状态变量。单输出信号y[n]的输出方程为：这里C是一个的1×k行向量（5-14）（5-15）22第二十二页，共五十七页，2022年，8月28日(5-16)23第二十三页，共五十七页，2022年，8月28日可算图在主对角线及其上方元素全部为零，具有这种属性的信号流图称为“可算图”：状态变量可以依次计算出来可算的充要条件：每一个闭合回路中的路径上至少含有一个延时单元，即z－1。24第二十四页，共五十七页，2022年，8月28日Maltab生成对于大量不同的模拟原型，可以轻易算出滤波器系数ak和bk。尽管生成的这些滤波器系数通常会用MATLAB子程序filter，但是使用下列三行简单的MATLAB程序，就可以逐个采样点地仿真IIR滤波器。25第二十五页，共五十七页，2022年，8月28日Maltab生成其中a和b在仿真循环外确定参数sreg代表长度为order+1的移位寄存器，其中order是滤波器阶数。可使用MATLAB滤波器综合子程序，如butter（巴特沃思滤波器）、cheby1（切比雪夫1型滤波器）及elliptic（椭圆滤波器），来计算滤波器系数向量a和b。必须小心以确保存放分子系数和分母系数的向量保持相同的长度，如果这两个向量长度不相等，较短的向量应该进行补零，使之与较长向量的长度相等。26第二十六页，共五十七页，2022年，8月28日例5-1来确定一个四阶巴特沃思滤波器的冲激响应。n=40; %numberofsamplesorder=4; %filterorder[b,a]=butter(order,0.1??); %prototype%Thefollowingprogramsegmentistheblockprocessingimplementation.in1=[1,zeros(1,n-1)]; %inputvectorout1=filter(b,a,in1); %outputvector%Thefollowingprogramsegmentisthesample-by-sampleimplementation.sreg=zeros(1,order+1); %initializeshiftregisterfork=1:nifk==1in=1; %impulseinputelsein=0;end27第二十七页，共五十七页，2022年，8月28日out=b(1)*in+sreg(1,1); %determineoutputsreg=in*b-out*a+sreg; %updateregistersreg=[sreg(1,2:(order+1)),0]; %shiftout2(k)=out; %createoutputvectorendsubplot(2,1,1)index=0:n-1;stem(index,out1)xlabel('SampleIndex')ylabel('BlockProcessing')subplot(2,1,2)stem(index,out2)xlabel('SampleIndex')ylabel('SerialProcessing')28第二十八页，共五十七页，2022年，8月28日Result29第二十九页，共五十七页，2022年，8月28日

5.3.2FIR滤波器的实现令ak＝0，可从图5-2直接得出FIR数字滤波器的实现策略。抽头延时线或横向延时线（TDL）结构，它们用带抽的延时线（级联式单元z－1）来实现，其中的抽头乘上了权值。30第三十页，共五十七页，2022年，8月28日

5.4IIR滤波器：综合方法与滤波器特性

引言：无限时间冲激响应（IIR）滤波器通常根据模拟原型（巴特沃思、切比雪夫、椭圆等）来进行设计，并采用递归结构实现。除射频端以外，DSP物理硬件的滤波器已经是数字形式的了，模拟原型为起点的数字滤波器综合的经典方法，就变得不太重要。根据模拟原型设计仍是基本和简单方法31第三十一页，共五十七页，2022年，8月28日

5.4.1冲激不变滤波器假设数字滤波器的单位冲激响应，除了幅度比例有差别外，等于模拟滤波器冲激响应的采样该滤波器称为模拟滤波器的冲激不变实现。滤波器的传递函数表示L－1拉普拉斯反变换，Z表示z变换，Ha(s)表示模拟原型的传递函数，T表示采样周期，同时还是幅度比例系数。（5-17）32第三十二页，共五十七页，2022年，8月28日一阶低通滤波器，其直流增益(f=0)为1。例5-2由式（5-17），首先对进行拉普拉斯反变换，由此得到如下冲激响应：（5-17）（5-18）u(t)为单位阶跃函数，表示所得的滤波器是因果的。进行z变换之前，对冲激响应在采样时刻进行采样。例子33第三十三页，共五十七页，2022年，8月28日z变换为进行求和并乘上采样周期得：的冲激不变等价(5-22)34第三十四页，共五十七页，2022年，8月28日5.4.2阶跃不变滤波器阶跃不变滤波器常用于仿真系统中阶跃不变滤波器:数字滤波器的单位阶跃响应等价于模拟滤波器阶跃响应的采样。模拟滤波器阶跃响应的频域表示为Ha(s)/s，其阶跃不变滤波器的传递函数如下：(5-23)35第三十五页，共五十七页，2022年，8月28日阶跃不变滤波器左边是传递函数为H(z)的数字滤波器对采样后的单位阶跃信号的响应，右边表示原型模拟滤波器单位阶跃响应采样的z变换。再次采用（5-18）式作为模拟原型。这个模拟原型的阶跃响应如下：使用部分分式展开可得(5-23)(5-24)(5-25)36第三十六页，共五十七页，2022年，8月28日取拉普拉斯反变换，并在t=nT处采样得再进行z变换得阶跃不变数字滤波器的传递函数为(5-26)(5-27)(5-28)OR(5-29)37第三十七页，共五十七页，2022年，8月28日Result将上式写成式（5-3）定义的标准形式有(5-30)38第三十八页，共五十七页，2022年，8月28日

5.4.3双线性z变换滤波器引言双线性z变换滤波器可能是仿真程序中使用最普遍的滤波器了采用双线性变换方法的数字滤波器综合很直观并且完全是代数形式的双线性变换不会表现出混叠现象。???39第三十九页，共五十七页，2022年，8月28日双线性z变换滤波器－

综合方法简单的代数变换，将模拟原型Ha（s)映射为一个数字滤波器，数字滤波器的传递函数定义为C为常数，因为Ha（s)通常定义为的多项式之比，（5-31）定义的代数映射刚好给出的合适形式，即的z-1多项式之比。具有理想的属性，why？exp(sdT)来替换双线性z变换中的变量z，sd是数字滤波器的复频域变量：（5-31）（5-32）40第四十页，共五十七页，2022年，8月28日双线性z变换滤波器－

综合方法cont1用j2πfa替换s,j2πfd替换sd（5-33）得到数字滤波器和模拟原型的关系调节C值使得对任意采样频率都有fa＝fd的方法称为预畸变（5-35）（5-34）或者41第四十一页，共五十七页，2022年，8月28日双线性z变换滤波器－

综合方法cont2fa

，fd显然是非线性关系，余切函数期望：模拟数字传递函数的形状保持不变尽可能线性，幅度相位不变x0,

近似线性变换要求fs是采样频率.满足的近似线性关系所需C值。（5-35）线性关系要求（5-36）如要求fa

＝fd（5-37）42第四十二页，共五十七页，2022年，8月28日例5-4为了说明双线性变换综合方法，我们再次使用（5-18）定义的简单一阶模拟原型。标准形式（5-18）（5-40）（5-39）（5-41）43第四十三页，共五十七页，2022年，8月28日例5-5这个例子将对前面三个例子中设计的三种滤波器的响应进行比较。假设采样频率为100Hz（Ts＝0.01s），参数为a=2π×10，从而使模拟滤波器的3dB频率为10Hz。冲激不变滤波器的传递函数为阶跃不变滤波器的传递函数为

（5-45）（5-43）（5-42）（5-44）44第四十四页，共五十七页，2022年，8月28日Howabout双线性z变换第一种如果不采用频率预畸变，那么模拟、数字在低频段非常吻合将C、a的值代入（5-48）（5-47）（5-46）45第四十五页，共五十七页，2022年，8月28日Howabout双线性z变换第二种设选择C值，使得频率响应在模拟滤波器的3dB频率处相匹配（5-51）（5-50）（5-49）（5-48）46第四十六页，共五十七页，2022年，8月28日三个数字滤波器比较%T=0.01;f=0:0.1:50;z=exp(-i*2*pi*f*T);%see(5.4)a0=0.239057;a1=0.239057;b1=0.521886;%bilinearinvariantnum=a0+a1*z;den=1-b1*z;ampx=abs(num./den);a0=0.628319;b1=0.533488; %impulseinvariantnum=a0;den=1-b1*z;ampy=abs(num./den);a0=1.0;a1=0.533488;b1=0.533488; %stepinvriant num=(a0-a1)*z;den=1-b1*z;ampz=abs(num./den);plot(f,ampx,f,ampy,f,ampz)xlabel('Frequency-Hz')ylabel('AmplitudeResponse')47第四十七页，共五十七页，2022年，8月28日Result单位增益零点48第四十八页，共五十七页，2022年，8月28日例5－6前面的例子为非常简单的一阶原型，ak,bk可以不用计算机确定考虑一个五阶椭圆滤波器，具有1dB通带波纹（PassbandRipple）和一个最小20dB的阻带衰减。也可以用前面的方法做，分析方法比较复杂也可以采用Matlab的综合工具49第四十九页，共五十七页，2022年，8月28日fs=100; %setsamplingfrequencyfc=20; %setcuttofffrequencyf=0:0.1:50; %definefrequencyvector[b,a]=ellip(5,1,20,2*pi*fc,'s'); %synthesizeellipticfilterh=freqs(b,a,2*pi*f); %amp.resp.ofanalogfilter[bz1,az1]=impinvar(b,a,fs);%impulseinvariantdigitalfilterh1=freqz(bz1,az1,f,fs); %amplituderesponseofabove[bz2,az2]=bilinear(b,a,fs); %bilinearzfilter(notprewarped)h2=freqz(bz2,az2,f,fs); %amplituderesponseofabove[bz3,az3]=bilinear(b,a,fs,fc); %bilinearzfilter(prewarped)h3=freqz(bz3,az3,f,fs); %amplituderesponseofabovesubplot(211)