通信系统的计算机模拟第七讲.ppt

1,通信系统的计算机模拟,第七讲,2,Announcement,Submit your proposal today,3,第六讲 回顾,IIR and FIR Filter IRR滤波器的实现 产生转置直接型滤波器的规则 矩阵形式 Maltab生成 FIR滤波器的实现 IIR滤波器：综合方法与滤波器特性 冲激不变滤波器 阶跃不变滤波器 双线性z变换 数字滤波器的计算机辅助设计,4,5.5 FIR滤波器：综合技术与滤波器特性,如果一个滤波器的冲激响应是有限长的，or 冲激响应截短为有限长的，则离散时域的滤波器输出信号为：,采用标准DSP表示法,（5.60）,（5.61）,FIR 模拟滤波器卷积积分的离散形式 （5-60）所定义的时域卷积运算可以用MATLAB函数filter进行仿真,（5.59）,5,FIR滤波器的吸引力,1通信中有的滤波器不能用拉普拉斯变换表示，即不能直接用IIR，如平方根升余弦（SQRC）脉冲整形滤波器和Jakes多普勒滤波器（Jakes Doppler filter）。FIR设计方法则可以轻易地实现这些滤波器的仿真。 2在许多仿真应用中，滤波器数据可能是以测量的频率响应或脉冲响应的形式经验性地给出， FIR仿真容易。（虽然有些方法可以拟合成AMMA自回归滑动平均，但频率响应数据不平滑时误差大） 3采用FIR滤波器，可以给定任意幅度响应和相位响应。 4FIR滤波器没有反馈，因此总是稳定的。,6,FIR仿真模型的不足,计算效率不如IIR滤波器实现的高。 （5-60）定义的卷积运算，产生一个输出采样点就需要N次（复数）乘法和加法运算，比如说N1024 点的话，FIR算法运行起来比典型的IIR滤波器算法 慢得多。 IIR的复杂度由滤波器的阶数决定的，而不同于FIR滤波器，由冲激响应的长度决定。,（5-60）,7,克服FIR仿真模型不足的措施,（5-60）,通过使用DFT/FFT算子来进行卷积运算，可以改进FIR模型的计算效率。注意： 1DFT/FFT运算具有固有的周期性，会产生圆周或周期卷积输出，为了线性卷积输出，输入和冲激响应进行填零补齐（Zero-padded）。2的整数次幂基-2FFT 2如果输入序列很长，通常把输入序列分成几个较小的互不重叠的子块，并对每个输入子块分别进行卷积，以得到对应的输出子块，然后将这些子块进行适当的交叠再相加。就得出总的输出。（matlab function fftfilt) 3基于DFT/FFT的卷积是一种分块处理运算。只有在输入采样点的数量累积到足够进行DFT/FFT运算之后，才可以产生输出子块的第一个采样点。不能用于仿真反馈循环中的滤波器。,8,时域卷积与基于DFT运算的卷积的计算效率,时域卷积与基于DFT运算的卷积的计算效率之比值跟log2N/N成比例。 如果冲激响应长度少于128个采样点的话，时域卷积与基于DFT/FFT的卷积在计算量上没有太大的区别。 如果对冲激响应的长度进行截短，FIR滤波器方法计算效率可以得到改进。截短等价于使用一个长度为N个采样的窗函数wk与冲激响应hk(k=1,2,.)相乘。 N值的选择？该窗函数内至少包含有总能量的98%。,9,窗函数,时域的截短（加窗）等价于频域的卷积,理想窗函数的频域表达式应该是 （f）。 频域理想时域不理想 窗函数的选择因此也是一种折中，所得窗口同时也引入了失真。 1频域的“窄主瓣”包含了大多数的能量。 2小旁瓣。 短形窗、汉明窗、Kaiser窗,最简单的矩形函数,（3-63）,10,矩形窗,11,窗函数,Hamming window,12,窗函数,Kaiser窗,I0 一类零阶修正贝塞尔函数， 确定窗函数形状的任意实数, N 给出窗长 (N + 1 点).,13,5.5.1 根据幅度响应进行的FIR设计,依据：滤波器的频率响应（包括幅度和相位）和单位冲激响应是一个傅里叶变换对。 推导单位冲激响应的一般做法：给定期望的幅度响应A(f)，再计算其傅里叶反变换 通常目标幅度响应为实偶函数，单位冲激响应也是实偶函数型的。因为冲激响应是偶函数，是非因果的. 如何变成时域上可实现的系统？ 必须对所得的冲激响应进行截短使其为有限长 问题及解决办法？截短冲激响应时，必须小心以避免引入显著的误差，适当的窗函数可以减少误差的影响 还要在时间轴上进行平移，使之成为因果函数 在时间轴上平移冲激响应，相当于使滤波器获得一个大小等于该时间平移长度的群时延的线性相位响应，可以得到具有任意幅度响应和线性相位响应的滤波器。 想要一个同时具有目标幅度响应和相位响应的滤波器，就要给定一个复传递函数。假定截断误差不显著，那么除了定常群时延外，最终得到的滤波器可以同时满足目标滤波器对幅度响应和相位响应的要求。,14,根据幅度响应进行的FIR设计方法,数字滤波器的幅度响应，是频率f的连续函数，是以采样频率为周期的周期函数，可表示为傅里叶级数。 这些离散的傅里叶级数的系数，给出期望的数字滤波器的冲激响应。 期望频率响应的IFT or期望频率响应采样点的IFFT冲激响应 这些方法很基本，小心使用会得到有用的滤波器设计,15,根据幅度响应进行的FIR设计方法,由于数字滤波器的幅度响应具有以采样频率为周期的周期性，其幅度响应可以展开为傅里叶级数，形式如下：,（5-64）,hn表示傅里叶级数的系数，M表示冲激响应的长度，而T=1/fs。,（5-4）,（5-6）,16,进一步化简,设M=2L+1，利用这一替换并改变序数k=n-L，可得：,h1k=hk+L,（5-68）,（5-67）,（5-66）,（5-65）,或者,17,分析,h1k=hk+L这是由式（5-64）所定义的因果滤波器的冲激响应平移L个采样点后得到的。 （5-68）定义滤波器不是一个因果滤波器。 但是，它便于采用基本的傅里叶变换法进行设计，并且只要简单地对冲激响应进行合适的平移就可以变为因果的。 很显然，由hn和h1n所确定的滤波器的幅度响应是相同的，而传递函数只是相差一个由式（5-67）所定义的线性相位平移。 因此，设计过程基于式（5-68）。,（5-68）,18,设计过程,假设一个滤波器将具有给定幅度响应,在式（5-68）两边同时乘以,（5-69）,在一个仿真带宽上，即 的一个周期内，（5-69）两边同时积分得,（570） （571）,这里,19,求得基本设计方程,这是我们的基本设计方程,（572）,（573）,代入（5-70）,20,例5-8 理想数字低通滤波器,我们设计了一个数字低通滤波器，来近似一个带宽为 的理想数字滤波器，这里fN为奈奎斯特频率，fs为采样频率，是介于0和1之间的一个参数。因此，期望的幅度响应为,（574）,21,解：,从（573） （575）,（576）,积分,由于fsT1，式（5-76）可以改写为,（577）,注意，这里,22,Matlab 实现,L = 30; % 2L+1 total points lam = 0.3; % normalized cutoff frequency m = -L:1:L; % vector of points bp = sin(pi*lam*(m+eps)./(pi*(m+eps); % impulse response为了防止当时出现不确%定形式，在序数中加上了eps?。 stem(0:2*L,bp,.) % plot impulse response xlabel(Sample index) ylabel(Impulse response),Matlab产生的第一个输出信号是冲激响应，h300.3。冲激响应关于mL=30处的中心是偶对称的。,23,Matlab 实现cont 1,figure; a = 1; freqz(bp,a) % plot amp and phase response,Freqz产生FIR滤波器的幅度响应和相位响应，再用默认的画图模式画出 线性相位相应，因为冲击响应关于中心偶对称 线性相位平移使滤波器具有因果性 幅度响应的符号变换会在阻带范围内导致锯齿状的相位响应。,24,Matlab 实现 cont2.,figure; subplot(2,1,1) % new figure H w = freqz(bp,a); plot(w/pi,abs(H); grid; % unwindowed amp response xlabel(Frequency (normalized to the Nyquist frequency = fs/2) ylabel(|H(f)| (unwindowed),线性标尺重新画出幅度响应，以避免对数坐标产生的幅度压缩。 波纹截短效应 矩形窗这里的加窗等于等价于把理想频率响应与 进行卷积（的数值由窗函数的宽度决定）。,25,Matlab 实现 cont3,subplot(2,1,2) w = 0.54+0.46*cos(pi*m/L); % Hamming window wbp = bp.*w; % apply window H w = freqz(wbp,a); plot(w/pi,abs(H); grid; % windowed amp response xlabel(Frequency (normalized to the Nyquist frequency = fs/2) ylabel(|H(f)| (windowed),通过采用一个从wn=1到wn=0较平滑过渡的窗函数，来降低矩形窗所造成的通带波纹。 汉明窗的权值定义如下,（5-78）,26,Matlab 实现 cont3,让意：窗函数必须在时间上平移L个采样点，从而使w0成为hn的中心。：采用汉明窗之后，通带/阻带波纹受到抑制,27,例5-9,设计了一个数字低通滤波器，不仅具有模拟巴特沃思滤波器的幅度响应，而且具有线性相位响应。我们没有该滤波器的等价模拟原型。巴特沃思滤波器通过如下幅度响应定义：,（5-79）,fc-带宽，3dB频率点，n-滤波器阶数,28,解：,2）然后对频率响应的采样值进行FFT反变换来产生冲激响应。 N点FFT，在序号（N/2)+1到N的范围内，返回的向量中会出现负时间采样点，这时必须重新编写序号，才能得到合适的冲激响应。,1）对频率响应进行采样,29,Matlab实现,order = 30; fc =5; % set filter parameters fmax = 100; % set max frequency npts = 256; % set number of samples f = (0:(npts-1)*(fmax/(npts-1); % frequency vector nn = 2*npts; % size ifft H = zeros(1,nn); % initialize vector Ha = 1./(sqrt(1+(f/fc).order); % amplitude response H = Ha 0 fliplr(Ha(2:npts); % even amplitude response cimp_resp = ifft(H,nn); % complex impulse response imp_resp = real(cimp_resp); % take real part aa = imp_resp(1:npts); % time = 0 bb = imp_resp(npts+1):nn); % time 0 reimpulse = bb aa; % real and even imp. resp. plot(reimpulse) % plot result,30,结果,接近理想sinc函数,31,Order=1,可以加窗忽略外面的值,32,5.5.2 根据冲激响应进行的设计,上例通过幅度响应和相位响应特性的指标来综合数字滤波器，是频率采样的一个应用。 在许多场合，已知滤波器冲激响应的解析表达式。 对冲激响应进行采样，就可以实现FIR滤波器的设计。 用于信号传输中成形脉冲的滤波器尤其如此 数据可以看作是在时间上以一定的脉冲周期分隔开的冲激序列。,（581）,dk1 or 1 in binary case 对应二进制的1 和 0 它通过一个具有冲激响应为p(t)的脉冲成形滤波器，得到：,（582）,33,例5-10 脉冲成形滤波,用于数据传输的脉冲波形p(t)，通常应选择具有奈奎斯特零符号间干扰（ISI）属性的脉冲波形。升余弦脉冲即是一种具有零ISI的脉冲波形，其中频域表达如下：,（583）,这里T是脉冲周期或者称为符号周期。,34,解：,对其进行傅里叶反变换，得到脉冲波形如下,（584）,非因果的脉冲 解决办法：通常将这个脉冲延时符号周期的整数倍，也就是mT，并将其长度截短到2mT，m的取值应该在便捷和精度要求之间进行折中。 然后在每个符号周期内进行点采样，使得TkTs，这里为采样周期。 用t-td=t-mT来进行替换t，令t=nTs ，及T=kTs得,（585）,35,求得pn,代入式（5-84），得到表示数字滤波器冲激响应的采样序列如下,（5-86）,令 为滤波器的输入，滤波器具有冲激响应,（5-87）,为了说明升余弦响应,36,Matlab实现,k = 10; % samples per symbol m = 4; % delay beta = 0.32; % bandwidth factor h = rcos(k,m,beta); % impulse response in = zeros(1,101); in(11) = 1; % input out = conv(in,h); % output t = 0:0.1:10; % time vector for plot stem(t,out(1:101),.) % plot output grid xlabel(Time) ylabel(Amplitude),37,rcos,function h=rcos(k,m,beta) % k - samples per symbol % m - delay is mT % beta - bandwidth factor beta = beta; n = 0:2*m*k; z = (n/k)-m+eps; t1 = cos(beta*pi*z); t2 = sin(pi*z)./(pi*z); t3 = 1-4*beta*beta*z.*z; h = t2.*t1./(t3);,38,对=0.32,k=10个采样点每符号，m=4个符号所得结果 脉冲周期整数倍处，pn=0,39,传递函数P(f)通常通过两个具有传递函数为 的滤波器的级联而实现。一个是发射器的一部分，另一个则是接收器的一部分。,（5-88）,使用式（5-85）定义的替换可实现延时和采样运算，所得结果为,（5-89）,平方根升余弦滤波器,40,冲激响应为Psqrcn的滤波器对输入信号n-1的响应,41,Matlab实现,T = 1; % symbol time k = 10; % samples per symbol m = 4; % delay beta = 0.32; % bandwidth factor h = sqrc(T,k,m,beta); % impulse response in = zeros(1,101); in(11) = 1; % input out = conv(in,h); % output h10 hn-1) out1 = conv(out,h); % conv of hn-10 and hn t = 2:0.1:8; % time vector for plot subplot(2,1,1),42,cont,stem(t,out(21:81),.) % plot hn-10 grid xlabel(Time) ylabel(hn-1) subplot(2,1,2) t = 6:0.1:12; % time vector for plot stem(t,out1(61:121),.) % plot conv of hn-10 and hn grid xlabel(Time) ylabel(conv(hn-1,hn),43,=0.32,k=10个采样点每符号，m=4符号,44,5.5.3 FIR滤波器仿真模型的实现,通信系统仿真中，FIR模型起着着键的作用 如果已知传递函数或者冲激响应，进行采样时的设计方法 FIR仿真模型的另外一个重要应用是在设计过程的后期阶段。滤波器也许已经设计并制造出来，这样就可以利用已经测量到的频率响应。 FIR滤波器仿真模型特别适用于频率响应以经验形式给定的滤波器仿真。,45,仿真步骤,前提：假定频率响应是任意的 表格形式经验性地给定（网络分析仪测量值获得） 第一步合理选择两个重要参数，即采样率和截短滤波器的持续时间. 准则： 采样率fs通常应满足fs16B到32B，这里B是滤波器带宽。 时间分辨率，采样间隔为 频率分辨率 理想的值处于B/64及B/32之间 采样点/符号数量必须为整数，通常是2的幂次方，最 小值为8。 冲激响应的持续时间通常选择8到16个符号。?,46,通常选取N=1024，N是滤波器抽头数量或者冲激响应的采样点数 确定了重要参数fs和N（从而也确定Ts和f）要对给定的频率响应数据进行预处理并实现仿真模型。 实现仿真模型的第一步涉及到对频率响应数据的预处理。 这一步骤包括：进行带通到低通的转换（如果这个滤波器是带通滤波器的话） 对频率响应数据进行重新采样，以及通过对群时延数据的数字积分将群时延转换为相位。 带通到低通的转换只要通过重新标记频率轴flffc就可以实现。 上述处理可导致低通等价滤波器不具有关于f0的共轭对称性。 低通滤波器的冲激响应可能取复数值 例子 残留边带滤波器低通转换的时候，就会产生这种影响。,47,第二步（重采样），诸如当网络分析仪测量结果所给出的频率点，跟要用于实现FIR滤波器的频率点不同时，这是必需的。 例如，网络分析仪收集到频率响应的100个采样点，而FIR滤波器的实现是基于N=1024个采样点的，对频率响应进行内插。简单的线性内插就足够了 最后一步群时延响应到相位响应的转换，通过在频域上对群时延进行数字积分，群时延数据就可以转换为相位响应 一旦得到低通等价滤波器的频率响应，使MATLAB实现FIR滤波器模型就包括以下步骤： 1将频率响应数据扩展到到之间fs/2,fs/2 ，并获得传递函数H(f)的采样值， ，这里 to N/21 ，f是频率分辨率。,2将H（f）的负频率部分平移到N/21与N之间，于是频率响应就包含在 H（kf） 之中，k1,2,3,.,N。,48,Cont.,3进行FFT反变换，就可以获得冲激响应时域采样率为 。 4如果必要的话，对冲激响应进行加窗（确保冲激响应作了合适旋转，且位于窗函数中心）。对加窗后的冲激响应进行归一化以便获得单位能量。 得到截短的冲激响应之后，使用a=1,b=h 的Matlab function filter来仿真滤波器，获得滤波器的输出信号。这里使用叠加法分块实现FIR滤波器，分块大小为N。,49,例511,卫星通信中使用的125MHz的带通信道滤波器。数据的测量特性由幅度响应和时延特性给出。采用前面讨论的方法进行FIR滤波器设计,50,51,52,Source file,fscale=1; fshift=0.0; dscale = 1000; % scaling parameters c5_Filter_Data; % load data Freq_Resp = data; fs = 900; filtsize = 512; ts = 1/fs; himp time = FIR_Filter_AMP_Delay(Freq_Resp,fs,filtsize,fscale,fshift,dscale); nw = 256; window1 = hamming(nw); window = zeros(filtsize,1); % Make sure the window is centered properly. wstart = (filtsize/2)-(nw/2); wend = (filtsize/2)+(nw/2)-1; window(wstart:wend) = window1; impw = himp.*window; figure; subplot(1,2,1); plot(abs(himp); grid; xlabel(Time Sample Index); ylabel(Filter Impulse Response); subplot(1,2,2); plot(abs(impw); grid; xlabel(Time Sample Index); ylabel(Windowed Filter Impulse Response); logpsd,freq,ptotal,pmax = log_psd(himp,filtsize,ts); lo