课 程 设 计课程设计名称: 数字信号处理 数字信号处理 专业课程设计任务书题 目 用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器主要内容 用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器,要求通带边界频率为400Hz,500Hz,阻带边界频率分别为350Hz,550Hz,通带最大衰减1dB,阻带最小衰减40dB,抽样频率为2000Hz,用MATLAB画出幅频特性,画出并分析滤波器传输函数的零极点;信号 经过该滤波器,其中 450Hz, 600Hz,滤波器的输出 是什么?用Matlab验证你的结论并给出 的图形。任务要求 1、掌握用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器的原理和设计方法。2、求出所设计滤波器的Z变换。3、用MATLAB画出幅频特性图。4、验证所设计的滤波器。参考文献 1、程佩青著,《数字信号处理教程》,清华大学出版社,20012、Sanjit K. Mitra著,孙洪,余翔宇译,《数字信号处理实验指导书(MATLAB版)》,电子工业出版社,2005年1月3、郭仕剑等,《MATLAB 7.x数字信号处理》,人民邮电出版社,2006年4、胡广书,《数字信号处理 理论算法与实现》,清华大学出版社,2003年审查意见 指导教师签字:说明:本表由指导教师填写,由教研室主任审核后下达给选题学生,装订在设计(论文)首页填 表 说 明1.“课题性质”一栏:A.工程设计;B.工程技术研究;C.软件工程(如CAI课题等);D.文献型综述;E.其它。2.“课题来源”一栏:A.自然科学基金与部、省、市级以上科研课题;B.企、事业单位委托课题;C.校、院(系、部)级基金课题;D.自拟课题。1 需求分析本课程设计要求用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器,给定的设计参数为通带边界频率为400Hz,500Hz,阻带边界频率分别为350Hz,550Hz,通带最大衰减1dB,阻带最小衰减40dB,抽样频率为2000Hz。要求采用双线性变换法,使得s平面与z平面是单值的一一对应关系,不存在频谱混淆的问题。由于数字频域和模拟频域的频率不是线性关系,这种非线性关系使得通带截止频率、过渡带的边缘频率的相对位置都发生了非线性畸变。设计出巴特沃兹型的带通滤波器之后,让信号 经过该滤波器,其中 450Hz, 600Hz,分析滤波器的输出 是什么,用Matlab验证结论并给出 的图形。2 概要设计1>用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点,可以采用非线性频率压缩方法,将整个频率轴上的频率范围压缩到-π/T~π/T之间,再用z=esT转换到Z平面上。也就是说,第一步先将整个S平面压缩映射到S1平面的-π/T~π/T一条横带里;第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系,消除了多值变换性,也就消除了频谱混叠现象,映射关系如图1-3所示。图1-3双线性变换的映射关系为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上,可以通过以下的正切变换实现(1-5)式中,T仍是采样间隔。当Ω1由-π/T经过0变化到π/T时,Ω由-∞经过0变化到+∞,也即映射了整个jΩ轴。将式(1-5)写成将此关系解析延拓到整个S平面和S1平面,令jΩ=s,jΩ1=s1,则得再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为:(1-6)(1-7)式(1-6)与式(1-7)是S平面与Z平面之间的单值映射关系,这种变换都是两个线性函数之比,因此称为双线性变换式(1-5)与式(1-6)的双线性变换符合映射变换应满足的两点要求。首先,把z=ejω,可得(1-8)即S平面的虚轴映射到Z平面的单位圆。其次,将s=σ+jΩ代入式(1-8),得因此由此看出,当σ<0时,|z|<1;当σ>0时,|z|>1。也就是说,S平面的左半平面映射到Z平面的单位圆内,S平面的右半平面映射到Z平面的单位圆外,S平面的虚轴映射到Z平面的单位圆上。因此,稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。2>双线性变换法优缺点双线性变换法与脉冲响应不变法相比,其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单值地对应于Z平面单位圆一周,即频率轴是单值变换关系。这个关系如式(1-8)所示,重写如下:上式表明,S平面上Ω与Z平面的ω成非线性的正切关系,如图7-7所示。由图7-7看出,在零频率附近,模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系;但当Ω进一步增加时,ω增长得越来越慢,最后当Ω→∞时,ω终止在折叠频率ω=π处,因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象,从而消除了频率混叠现象。图1-4双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的,如式(1-8)及图1-4所示。由于这种频率之间的非线性变换关系,就产生了新的问题。首先,一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器,不再保持原有的线性相位了;其次,这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的,即某一频率段的幅频响应近似等于某一常数(这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性),不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变,如图1-5所示。图1-5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器,双线性变换后,仍得到幅频特性为分段常数的滤波器,但是各个分段边缘的临界频率点产生了畸变,这种频率的畸变,可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变,然后经变换后正好映射到所需要的数字频率上。根据以上IIR数字滤波器设计方法,下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器,通带截止频率fp1=500HZ,fp2=400HZ;阻带截止频率fs1=350HZ,fs2=550HZ;通带最大衰减αp=1dB;阻带最小衰减αs=40dB,抽样频率为2000HZ。I.设计步骤:(1)根据任务,确定性能指标:通带截止频率fp1=500HZ,fp2=400HZ;阻带截止频率fs1=350HZ,fs2=550HZ;通带最大衰减αp=1dB;阻带最小衰减αs=40dB;抽样频率为2000HZ。 (2)用Ω=(2/T)tan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变,得到带通模拟滤波器H(s)的边界频率主要是通带截止频率fp1,fp2;阻带截止频率fs1,fs2的转换。为了计算简便,对双线性变换法一般T=2s通带截止频率fc1=(2/T)*tan(2*pi*fp1/2)fc2=(2/T)*tan(2*pi*fp2/2)阻带截止频率fr1=(2/T)*tan(2*pi*fs1/2)fr2=(2/T)*tan(2*pi*fs2/2)阻带最小衰减αs=1dB和通带最大衰减αp=40dB;(3)运用低通到带通频率变换公式λ=(((f^2)-(f0^2))/ f)将模拟带通滤波器指标转换为模拟低通滤波器指标。B=fc2-fc1normwr1=(((fc1*fc2)-(fc1^2))/fc1) normwr2=(((fc1*fc2)-(fc2^2))/fc2)normwc1=(((fc1*fc2)-(fr1^2))/fr1)normwc2=(((fc1*fc2)-(fr1^2))/fr1)模拟低通滤波器指标:normfc,normfr,αp,αs(4)设计模拟低通原型滤波器。借助巴特沃斯(Butterworth)滤波器,用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数H(s);(5)调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。(6)利用双线性变换法将模拟带通滤波器H(s)转换成数字带通滤波器H(z).II.程序流程框图:开始↓读入数字滤波器技术指标↓将指标转换成归一化模拟低通滤波器的指标↓设计归一化的模拟低通滤波器阶数N和1db截止频率↓模拟域频率变换,将G(P)变换成模拟带通滤波器H(s)↓用双线性变换法将H(s)转换成数字带通滤波器H(z)↓输入信号后显示相关结果↓结束3 运行环境PC机 MATLAB4 开发工具和编程语言MATLAB语言5 详细设计clc;clear all;%设计滤波器%所需设计的带通滤波器的参数设置Fres=2000;Ts=1/Fres;Omgap1=500Omgap2=400Omgas1=350Omgas2=550%进行双线性变换,使得s平面与z平面是单值的一一对应关系Omgap=(Omgap1-Omgap2)*TsOmgap3=tan(Omgap1*Ts*2*pi/2)Omgap4=tan(Omgap2*Ts*2*pi/2)Omgas3=tan(Omgas1*Ts*2*pi/2)Omgas4=tan(Omgas2*Ts*2*pi/2)%对参数归一化ap1=Omgap3/Omgapap2=Omgap4/Omgapas1=Omgas3/Omgapas2=Omgas4/Omgap%将模拟带通滤波器指标转化成模拟低通滤波器指标I1=(ap1*ap2-as1*as1)/as1I2=(ap1*ap2-as2*as2)/as2I3=(ap1*ap2-ap1*ap1)/ap1I4=(ap1*ap2-ap2*ap2)/ap2I5=min(I1,-I2)alfpp=1;alfps=40;%设计巴特沃兹型低通滤波器HL(s)[N,Wn]=buttord(I4,I5,alfpp,alfps,'s')[Num,Den]=butter(N,Wn,'s')%将设计的模拟低通滤波器传递函数HL(s)转化成模拟带通滤波器H(s)[Num2,Den2]=lp2bp(Num,Den,sqrt(Omgap3*Omgap4),Omgap);%将设计的模拟带通滤波器H(s)转化成数字带通滤波器H(Z)[Num3,Den3]=bilinear(Num2,Den2,0.5);%画出H(Z)幅频特性曲线w=0:pi/255:pi;h=freqz(Num3,Den3,w);H=20*log10(abs(h));plot(w,H);%设计信号函数f1=450;f2=600;t=0:1/2000:40/f2;f=sin(2*pi*f1*t)+sin(2*pi*f2*t);plot(f);%信号通过带通滤波器g=invfreqz(h,w,40,50)g1=conv(g,f)plot(g1)6 调试分析在设计滤波器的时候,计算的是幅频特性。而相频特性没有进行频谱分析。而信号和信号通过滤波器的图形分析都为时与分析,没有进行频域分析,是不知道该用哪个函数对其进行傅里叶变换。试过freqz,fft等函数,但是都没有出来结果。说明调用格式不正确。本次试验用的是巴特沃兹型低通滤波器设计的方法,还可以利用切比雪夫型低通滤波器设计,或者椭圆型等。在本次试验中,我们只用了两种频率简单的叠加作为输入信号,但实际信号是多种频率的组合,可以再增加一些频率。现实中还应该有噪声的影响,本实验中没有考虑。可以再加上噪声信号。7 测试结果Fres = 2000Ts = 5.0000e-004Omgap1 = 500Omgap2 = 400Omgas1 = 350Omgas2 = 550Omgap = 0.0500Omgap3 = 1.0000Omgap4 = 0.7265Omgas3 =0.6128Omgas4 = 1.1708ap1 = 20.0000ap2 = 14.5309as1 = 12.2560as2 = 23.4170I1 = 11.4562I2 = -11.0065I3 = -5.4691I4 = 5.4691I5 = 11.0065alfpp = 1alfps = 40N = 8Wn = 6.1894Num = 1.0e+006 * Columns 1 through 8 0 0 0 0 0 0 0 0 Column 9 2.1538Den = 1.0e+006 * Columns 1 through 8 0.0000 0.0000 0.0005 0.0052 0.0377 0.1984 0.7386 1.7837 Column 9 2.1538Num2 = 1.0e-004 * Columns 1 through 8 0.8413 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Column 9 -0.0000Den2 = Columns 1 through 8 1.0000 1.5863 7.0705 8.7151 20.5005 19.9985 32.1353 24.8474 Columns 9 through 16 29.9185 18.0527 16.9631 7.6698 5.7123 1.7643 1.0400 0.1695 Column 17 0.0776Num3 = 1.0e-004 * Columns 1 through 8 0.0043 0.0000 -0.0341 0.0000 0.1194 0.0000 -0.2389 0.0000 Columns 9 through 16 0.2986 0.0000 -0.2389 0.0000 0.1194 0.0000 -0.0341 0.0000 Column 17 0.0043Den3 = Columns 1 through 8 1.0000 -2.2463 8.3946 -13.4519 27.6744 -33.6694 48.3386 -45.7295 Columns 9 through 16 49.5687 -36.4140 30.6576 -16.9904 11.1207 -4.2944 2.1332 -0.4524 Column 17 0.1603h = Columns 1 through 5 -0.0000 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i …… Columns 251 through 256 -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000iH = Columns 1 through 8 -279.2737 -279.2732 -279.2723 -279.2813 -279.3855 -280.0020 -283.0382 -292.4507…… Columns 249 through 256 -281.0032 -280.3352 -280.1410 -280.0979 -280.0943 -280.0973 -280.0996 -280.1004f1 = 450f2 = 600t = Columns 1 through 8 0 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 …… Columns 129 through 134 0.0640 0.0645 0.0650 0.0655 0.0660 0.0665f = Columns 1 through 8 0 1.9387 -0.2788 -1.4788 0.3633 0.7071 -0.1420 0.1338…… Columns 129 through 134 -0.3633 -0.7946 1.0000 1.1075 -1.5388 -1.0418g = Columns 1 through 8 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0001 -0.0002…… Columns 33 through 41 -0.0002 0.0000 0.0001 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000g1 = Columns 1 through 8 0 0.0000 0.0000 -0.0000 -0.0000 0.0001 0.0001 -0.0003…… Columns 169 through 174 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000综合上述分析,数据基本正确。上图为带通滤波器幅频响应 、上图为传输函数零极点分析 上图为输入信号时域上图为输出波形的时域分析参考文献 [1] 吴大正 《信号与线性系统分析》第四版 高等教育出版社[2] 郑君里 《信号与系统》第二版 高等教育出版社[3] Sanjit K. Mitra 《数字信号处理—基于计算机的方法》第三版 清华大学出版社[4] 余成波 《数字信号处理及MATLAB实现》清华大学出版社[5] 周利清 《数字信号处理基础》 北京邮电大学出版社[6] 美国莱昂 《数字信号处理》英文第二版 机械工业出版社心得体会................0.0
摘 要 FIR数字滤波器是数字信号处理的经典方法,其设计方法有多种,用DSP芯片对FIR滤波器进行设计时可以先在MATLAB上对FIR数字滤波器进行仿真,所产生的滤波器系数可以直接倒入到DSP中进行编程,在编程时可以采用DSP独特的循环缓冲算法对FIR数字滤波器进行设计,这样可以大大减少设计的复杂度,使滤波器的设计快捷、简单。关键词 FIR;DSP;循环缓冲算法1 引言在信号处理中,滤波占有十分重要的地位。数字滤波是数字信号处理的基本方法。数字滤波与模拟滤波相比有很多优点,它除了可避免模拟滤波器固有的电压漂移、温度漂移和噪声等问题外,还能满足滤波器对幅度和相位的严格要求。低通有限冲激响应滤波器(低通FIR滤波器)有其独特的优点,因为FIR系统只有零点,因此,系统总是稳定的,而且容易实现线性相位和允许实现多通道滤波器。2 FIR滤波器的基本结构及设计方法2.1 FIR滤波器的基本结构设a i(i=0,1,2,…,N一1)为滤波器的冲激响应,输入信号为 x(n),则FIR滤波器的输入输出关系为: FIR滤波器的结构如图1所示:图12.2 FIR滤波器的设计方法 (1) 窗函数设计法 从时域出发,把理想的无限长的hd(n)用一定形状的窗函数截取成有限长的h(n),以此h(n)来逼近hd(n),从而使所得到的频率响应H(ejω)与所要求的理想频率响应Hd(ejω) 相接近。优点是简单、实用,缺点是截止频率不易控制。 (2) 频率抽样设计法从频域出发, 把给定的理想频率响应Hd(ejω)以等间隔抽样,所得到的H(k)作逆离散傅氏变换,从而求得h(k),并用与之相对应的频率响应H(ejω)去逼近理想频率响应Hd(ejω)。优点是直接在频域进行设计,便于优化,缺点是截止频率不能自由取值。(3) 等波纹逼近计算机辅助设计法前面两种方法虽然在频率取样点上的误差非常小,但在非取样点处的误差沿频率轴不是均匀分布的,而且截止频率的选择还受到了不必要的限制。因此又由切比雪夫理论提出了等波纹逼近计算机辅助设计法。它不但能准确地指定通带和阻带的边缘,而且还在一定意义上实现对所期望的频率响应实行最佳逼近。3 循环缓冲算法对于N级的FIR滤波器,在数据存储器中开辟一个称之为滑窗的N个单元的缓冲区,滑窗中存放最新的N个输入样本。每次输入新的样本时,一新样本改写滑窗中的最老的数据,而滑窗中的其他数据不需要移动。利用片内BK(循环缓冲区长度)寄存器对滑窗进行间接寻址,环缓冲区地址首位相邻。下面,以N=5的FIR滤波器循环缓冲区为例,说明循环缓冲区中数据是如何寻址的。5级循环缓冲区的结构如图所示,顶部为低地址。……由上可见,虽然循环缓冲区中新老数据不很直接明了,但是利用循环缓冲区实现Z-1的优点还是很明显的:它不需要数据移动,不存在一个极其周期中要求能进行一次读和一次写的数据存储器,因而可以将循环缓冲区定位在数据存储器的任何位置(线性缓冲区要求定位在DARAM中)。实现循环缓冲区间接寻址的关键问题是:如何使N个循环缓冲区单元首位相邻?要做到这一点,必须利用BK(循环缓冲器长度)器存器实现按模间接寻址。可用的指令有:… *ARx+% ;增量、按模修正ARx:addr=ARx,ARx=circ(ARx+1)… *ARx-% ;减量、按模修正ARx:addr=ARx,ARx=circ(ARx-1)… *ARx+0% ;增AR0、按模修正ARx:addr=ARx,ARx=circ(ARx+AR0)… *ARx-0% ;减AR0、按模修正ARx:addr=ARx,ARx=circ(ARx-AR0)… *+ARx(lk)% ;加(lk)、按模修正ARx:addr=circ(ARx+lk),ARx=circ(ARx+AR0)其中符号“circ”就是按照BK(循环缓冲器长度)器存器中的值(如FIR滤波其中的N值),对(ARx+1)、(ARx-1)、(ARx+AR0)、(ARx-AR0)或(ARx+lk)值取模。这样就能保证循环缓冲区的指针ARx始终指向循环缓冲区,实现循环缓冲区顶部和底部单元相邻。循环寻址的算法可归纳为:if 0 index + step < BK: index = index + stepelse if index + step BK: index = index + step – BKelse if index + step < BK: index = index + step + BK上述算法中,index是存放在辅助寄存器中的地址指针,step为步长(亦即变址值。步长可正可负,其绝对值晓予或等于循环缓冲区长度BK)。依据以上循环寻址算法,就可以实现循环缓冲区首位单元相邻了。 为了使循环缓冲区正常进行,除了用循环缓冲区长度寄存器(BK)来规定循环缓冲区的大小外,循环缓冲区的起始地址的k个最低有效位必须为0。K值满足2k>N,N微循环缓冲区的长度。4 FIR滤波器在DSP上的实现对于系数对称的FIR滤波器,由于其具有线性相位特征,因此应用很广,特别实在对相位失真要求很高的场合,如调制解调器(MODEM)。例如:一个N=8的FIR滤波器,若a(n)=a(N-1-n),就是对称FIR滤波器,其输出方程为:y(n)= a0x(n)+ a1x(n-1)+ a 2x(n-2)+ a 3x(n-3)+ a 3x(n-4)+ a 2x(n-5)+ a1x(n-6)+ a0x(n-7)总共有8次乘法和7次加法,如果改写成: y(n)= a0 [x(n)+ x(n-7)]+ a1 [ x(n-1)+ x(n-6)]+ a 2 [ x(n-2)+ x(n-5)]+ a 3 [ x(n-3)+ x(n-4)]则变成4次乘法和7次加法。可见,乘法运算的次数减少了一半。这是对称FIR的又一个优点。对称FIR滤波器C54X实现的要点如下:(1)数据存储器中开辟两个循环缓冲算区:新循环缓冲区中存放新数据,旧循环缓冲区中存放老数据。循环缓冲区的长度为N/2。 (2)设置循环缓冲区指针:AR2指向新循环缓冲区中最新的数据,AR3指向旧循环缓冲区中最老的数据。 (3)在程序存储器中设置系数表。 (4)AR2+ AR3 AH(累加器A的高位),AR2-1AR2,AR3-1 AR3 (5)将累加器B清零,重复执行4次(i=0,1,2,3):AH*系数ai+B B,系数指针(PAR)加1。AR2+ AR3AH,AR2和AR3减1。 (6)保存和输出结果。 (7)修正数据指针,让AR2和AR3分别指向新循环缓冲区中最老的数据和旧循环缓冲区中最老的数据。 (8)用新循环缓冲区中最老的数据替代旧循环缓冲区中最老的数据,旧循环缓冲区指针减1。 (9)输入一个新的数据替代新循环缓冲区中最老的数据。 重复执行第(4)至(9)步。 在编程中要用到FIRS(系数对称有限冲击响应滤波器)指令,其操作步骤如下: FIR Xmem,Ymem,Pmem 执行 Pmad PAR 当(RC)≠0 (B)+(A(32-16))×(由PAR寻址Pmem)B ((Xmem)+(Ymem))<<16A (PAR)+1PAR (RC)-1RC FIRS指令在同一个及其周期内,通过C和D总线读2次数据存储器,同时通过P总线读一个系数 本文对FIR滤波器在DSP上的实现借助了MATLAB,其设计思路为:(1)MATLAB环境下产生滤波器系数和输入的数据,并仿真滤波器的滤波过程,可视化得到滤波器对动态输入数据的实时滤波效果;(2)将所得滤波器系数直接导入CCStudio中,再把滤波器的输入数据作为CCStudio设计的滤波起的输入测试数据存储在C54x数据空间中; (3)在CCStudio环境下结合FIR滤波的公式适用汇编语言设计FIR滤波程序,使用MATLAB产生的滤波器系数和输入测试数据进行计算,把输入数据和滤波结果借助CCStudio菜单中的View/Graph/Time/Frequency子菜单用图形方式显示出来(结果如图2);图2 (a)输入数据(Input)图2(b)滤波后的数据(Output) 将FIR滤波的入口数据地址改为外部I/O空间或McBSP口的读写数据地址,或数据空间内建缓冲地址;将FIR滤波的结果数据地址改为外部I/O空间或McBSP口的输出数据地址,或数据空间内建缓冲地址,则完成了基于C54xDSP的实时数据FIR滤波程序。参考文献:[1] 程佩青.数字信号处理教程[M].北京:清华大学出版社 1999年[2] 孙宗瀛,谢鸿林.TMS320C5xDSP原理设计与应用[M].北京:清华大学出版社.2002年[3] 陈亚勇等 编著.MATLAB信号处理详解[M].北京:人民邮电出版社.2001年[4] Texas Instruments.TMS320C54x Assembly Language Tools User’s Guide[5] Texas Instruments.TMS320C54x DSP Programmer’s Guide
85 浏览 3 回答
271 浏览 3 回答
146 浏览 1 回答
219 浏览 3 回答
292 浏览 5 回答
326 浏览 3 回答
271 浏览 4 回答
146 浏览 2 回答
138 浏览 2 回答
332 浏览 2 回答
349 浏览 2 回答
273 浏览 3 回答
263 浏览 1 回答
104 浏览 3 回答
269 浏览 2 回答