⑷ 最优线性相位FIR滤波器的设计步骤

（1）输入部分，包括滤波器性能要求及滤波器类型，前者指的是所需的频率响应的幅度函数，加权函数和滤波器单位抽样响应长度N,后者是要指出所需设计的是带通（包括低通、带通、高通、带阻等）滤波器或是微分器或是希尔伯特变换器；

（2）用公式表示逼近问题，也就是用式（3.3.5）来表示加权逼近误差;

（3）用瑞米兹多次交换算法，求逼近问题的解；

（4）计算滤波器的单位抽样响应。

这里步骤（1）表示设计者利用算法时必须具有的条件，步骤（2）在前面加权切比雪夫逼近中已讨论过，而步骤（3）正是瑞米兹交换算法求逼近问题的解。在整个设计程序中，瑞米兹交换算法是以子程序形式出现的。

更仔细地说，设计FIR数字滤波器的最优化算法包括：① 输入数据、滤波器性能要求及滤波器类型；② 根据滤波器的类型和单位抽样响应长度，确定逼近函数的个数r；③ 在从0到的频率区间，用密集的格点来表示离散频率，两格点的距离可表示成

将格点频率用有下标的频率来表示并赋以标称频率值；④ 用子程序EFF和WATE分别计算各格点频率上要求的函数值和加权函数值；⑤ 用公式逼近问题,即将式（3.3.4）转变为（3.3.4）式，将和分别变成和；⑥ 瑞米兹交换算法求逼近问题，这包括设定(r+1)个极值频率的初始假设值（猜想值）；⑦ 计算滤波器单位抽样响应；⑧ 给出仿真结果并把最佳误差和单位抽样响应打印出来。

下面简要介绍一下瑞米兹交换算法：

Remez 算法是由Parks 和McClellan 等人在1972年推导出来的。它是将FIR DF 的五个参数( N , ,,,) 中的N , ,和固定,而视(或)为变量的一种迭代方法。尽管Herrmann 等人在1971年也推导出来了将N 和,固定，而将,设为变量的另一种迭代方法, 但由于前一种算法最灵活且最有效。因而成为最优化设计的主要方法。该方法求解过程为：① 求解δ:首先在滤波器的通带和阻带内等间隔地取r + 1 个频率点( k = 0 ,1…r) 作为交错点的初始值, 然后利用式（3.3.7）求解出满足下式的δ值: 

， k=0,1,…,r

② 求解：利用求解出来的δ值和预先假定的r +1 个频率点求出(其中i = 0 ,1Λr - 1)的值,然后根据拉格朗日插值公式求出的最终表达式。③ 求解 :将求得的P(ω) 代入下式判定其是否满足不等式

=[-]

的要求。若满足要求，则说明已经获得了最优解； 若在某些频率点不满足要求，则需要将这些频率点作为新的极值点重新计算, 经过反复的迭代直到在所有的频率点上都满足不等式的要求。这时的值就是最终所求的值。这样便获得了最佳逼近。

                             滤波器类型

图3.3.1 最佳线性相位FIR滤波器设计算法框图

                                   Y

                  N

     Y

                N

图3.3.2 瑞米兹交换算法流图

4．FIR数字滤波器的MATLAB实现

4.1.MATLAB概述

MATLAB的历史  MATLAB的首创者是在数值线性代数领域颇有影响的Cleve Moler博士，他在讲授线性代数课程时，深感高级语言编程的诸多不便之处，于是萌生了开发新的软件平面，即为MATLAB(MATtix LABoratory,矩阵实验室)，软件采用了当时流行的EISPACK（基于特征值计算的软件包）和LINPACK（线性代数软件包）中的子程序，利用FORTRAN语言编写而成。现今的MATLAB已全部采用C语言编写，并使用户界面变得越来越好[1]。

由Cleve Moler博士等一批数学家和软件专家组建了Math Works软件公司，专门从事MATLAB的扩展和改进。自1982年推出第一个版本以来，1992年推出了具有划时代意义的MATLAB V4.0，1993年推出了可用于IBM PC及其兼容机上的微机版，特别是与Windows配合使用，使MATLAB的应用得到了前所未有的发展。1994年推出了成熟的4.2版本，并得到了广泛的重视和应用。MATLAB5.0版本在下列五个方面对以前版本进行了改进：

⑴ 增强了编程和应用开发工具；

⑵ 采用了新的数据类型、结构及语言特性；

⑶ 采用了更快、更好的图形和可视化技术；

⑷ 提供了更多数学和数据分析工具；

⑸ 对MATLAB应用工具箱和Simulink作了较大的增强。

MATLAB 5.1版本再次增强了系统能力，除了对5.0版本许多函数的改进和增强外，还提供了Stateflow 1.0,Mapping Toolbox 1.0,Real-Time Workshop 2.1等重要部件。

MATLAB 5.2再次对系统进行了改进。MATLAB 5.3版本已经成为一个功能比较强大、性能稳定的软件。

笔者此次设计采用的是目前最新的MATLAB 6.5版本，它的功能更加强大，使用更加方便。

MATLAB的特点  MATLAB之所以为广大读者所喜爱，是因为它具有其它语言所不具备的特点。

⑴ 在MATLAB中，以复数矩阵作为基本编程单元，使矩阵操作变得轻而易举。MATLAB中矩阵操作如同其它高级语言中的变量操作一样方便，而且矩阵无需定义即可采用，可随时改变矩阵的尺寸，这在其它高级语言中是很难办到的。

⑵ MATLAB语句书写简单，表达式的书写如同在稿纸中演算一样，与人们的手工运算一致，容易为人们所接受。

⑶ MATLAB语句功能强大，一条语句往往相当其它高级语言中的几十条，几百条甚至几千条语句。例如用MATLAB求解FFT问题时，仅需几条语句，而采用C语言实现时需要几十条，采用汇编语言实现则需要3000多条语句。

⑷ MATLAB系统具有丰富的图形功能。MATLAB系统本身是一个Windows下的具有良好用户界面的系统，而且提供了丰富的图形界面设计函数。如专门用于绘制二维曲线的plot函数，用于绘制三维曲线的plot3函数。在工具箱中，有些函数本身可以提供良好的图形功能，如step函数可计算指定系统的单位阶跃响应，并直接在屏幕窗口中绘制出系统的单位阶跃响应曲线。

⑸ MATLAB提供了许多面向应用问题求解的工具箱函数，从而大大方便了各个领域专家学者的使用。目前，MATLAB提供了20多个工具箱函数，如信号处理，图像处理，控制系统，非线性控制设计，鲁棒控制，系统辨识，最优化，神经网络，模糊系统和小波等，它们提供了各个领域应用求解的便利函数，使系统分析与设计变得更加简捷。

⑹ MATLAB的易扩展性是最重要的特性之一，也是MATLAB得以广泛应用的原因之一。MATLAB给用户提供了广阔的扩展空间，用户可以很容易编写出适用与自己和专业的M文件，供自己或同伴使用，这实际上就是扩展了MATLAB的系统功能[3]。

4.2．滤波器设计的MATLAB实现

MATLAB是一个功能强大的工具软件，其中Signal Processing Toolbox提供了2个配合使用的函数进行FIR滤波器的最优等波纹设计：remez和remezord。其调用格式有以下几个版本：

    ① remez:      

[h]=remez(N,f,m) 

[h]=remez(N,f,m,weights) 

[h]=remez(N,f,m,`ftype`) 

[h]=remez(N,f,m,weights,`ftype`) 

现分别说明如下：

（1）[h]=remez(N,f,m)设计一N阶（注意，滤波器长度是M=N+1）FIR数字滤波器，其频率响应由数组f和m给定。滤波器系数（或脉冲响应）在长度为M的数组h中得到。数组f包含以为单位的频率边缘频率，也即0.0≤f≤1.0，这些频率必须是以递增的次序从0.0开始，并以1.0终结。数组m包含有在f所给定频率上的期望幅度响应。f和m数组的长度必须相同而且必须是偶数。在每个频带内所用的加权函数都是等于1，这意味着在每一频带内的容度（）是相同的。

（2）[h]=remez(N,f,m,weights)除了数组weights给出了在每个频带内的加权函数外，其余都与上面情况是类似的。

（3）[h]=remez(N,f,m,ftype)与第一种情况类似，除了当ftype是字符串‘differentiator’或‘hilbert’时，它分别设计数字微分器和数字希尔伯特变换器。对于数字希尔伯特变换器，在数组f中的最低频率不应该是0，而最高频率也不应该是1。对数字微分器来说，m向量不给出在每个频带内的期望斜率而是期望幅度。

（4）[h]=remez(N,f,m,weights,ftype)除了数组weights给出每个频带的加权函数外，与上面情况类似。

正如在描述Parks－McClellan算法过程中所说明的，为了利用remez子程序必须首先利用（3.3.9）式估猜滤波器的阶。在得到的数组h中的滤波器系数之后，必须要校正最小阻带衰减并将它与给定的比较，然后增大（或降低）这个滤波器的阶。必须将这一过程重复直到得到所期望的为止。

② remezord

      [n,fo,mo,w]=remezord(f,m,dev)

      [n,fo,mo,w]=remezord(f,m,dev,Fs)

说明：

remezord函数可为remez选择滤波器的阶。给定频域中的性能指标，remezord可以产生近似满足指标的最小阶。

[n,fo,mo,w]=remezord(f,m,dev)可找出近似的阶n、归一化频带边界fo、频带内幅值mo及加权矢量w，使由remez函数构成的滤波器满足f,m及dev指定的性能要求。其中f和m用来指示频段及相应的幅值，这与firls函数类似，但应该注意一点，f中省去了1和0这两点，即

       length(f)=2*length(m)-2

dev用于指定各频段允许的偏差。

[n,fo,mo,w]=remezord(f,m,dev,Fs)可指定取样频率Fs，如果缺省Fs，则默认为2 Hz。

应该注意，有时remezord函数对阶估计不足，因此可能会使滤波器达不到指定的性能，这时应稍微增加阶次。

4.3．设计举例

例1．设计一17阶的Parks－McClellan带通滤波器，并画出期望幅频曲线和实际幅频曲线。

其程序如下：

f=[0 0.3 0.4 0.6 0.7 1];

m=[0 0 1 1 0 0];

b=remez(17,f,m);

[h,w]=freqz(b,1,512);

plot(f,m,w/pi,abs(h))

 这可得如下图所示的幅频曲线。

图4.3.1 理想和实际带通滤波器幅频曲线

例2．设计一最小阶的低通滤波器，通带截止频率为500Hz，阻带截止频率为600Hz(取样频率为2000Hz),通带波纹小于3dB，阻带波纹低于40dB。

其程序如下：

Rp=3; Rs=40;

Fs=2000; f=[500 600]

m=[1 0];

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1) 10^(-Rs/20)];

[n, fo, mo, w]=remezord(f, m, dev, Fs);

n=n+1;

b=remez(n, fo, mo, w);

[h, f]=freqz(b, 1, 1024, Fs);

plot(f, 20*log10(abs(h)))

               图4.3.2 低通滤波器特性（n=23）

4.4．设计任务

本次设计的任务是利用MATLAB设计出低通滤波器，高通滤波器，带通滤波器各一个。并且要求研究III型滤波器，由于其90度的相移，我又设计了长度为51的希尔伯特变换器一个。现给出各滤波器的指标要求，具体的原程序请参阅附录。

a．低通滤波器

通带截止频率＝0.2；阻带截止频率＝0.3；

通带容度＝0.25dB；   阻带容度＝50dB；

b．高通滤波器

通带截止频率＝0.75；阻带截止频率＝0.6；

通带容度＝0.5dB；     阻带容度＝50dB；

c．带通滤波器

通带截止频率＝0.35；＝0.65；

阻带截止频率＝0.2； ＝0.8；

通带容度＝1.0dB；阻带容度＝60dB。

程序的编写过程中，除了使用remez函数来计算滤波器的h(n)和长度M，还要用到freqz_m函数来计算滤波器的频率响应，freqz_m函数是freqz函数的修正形式，在调用前要先对其加以定义。

下面给出各滤波器时域和频域的波形图。

                         图4.4.1 低通滤波器

运行结果滤波器长度为47，阻带容度51.0857分贝。

图4.4.2 高通滤波器

运行结果滤波器长度为29，阻带容度50.2253分贝。

                      图4.4.3 希尔伯特变换器

运行结果长度51，处于频带之内。

图4.4.4 带通滤波器

运行结果滤波器长度29，阻带容度61.2843分贝。

    从运行的情况可以看出，滤波器的长度适中，通带和阻带的波纹均达到了设计指标要求。

5．总结

     FIR数字滤波器的三种设计方法可以说各有各的优势。窗函数法的优点是简单，有闭合形式的公式可循，因而很实用。其缺点是通带、阻带的截止频率不易控制。频率抽样法也可以对过渡带抽样进行优化设计，但它只是将过渡带的几个抽样点作为变量，在某一优化准则下，通过计算机进行迭代运算，以得到最优的结果，所以它只是接近最优化，而不是最优化。等波纹最佳逼近法虽然可以得到最优化的结果，但其使用的最大误差最小化准则和交错定理均不易理解，在求解逼近问题的解时，若采用人工计算，计算量将非常庞大，而且容易出错。但这些缺点均可由电脑来弥补，尤其用MATLAB实现比较方便。

MATLAB系统供了许多工具箱（Toolbox），借助于信号处理工具箱（signal processing）中的freqz_m，remez等函数，使得FIR数字滤波器的设计大为简化，每个程序都只有短短的几十行。从运行结果和实验波形可以看到设计出来低通、高通、带通滤波器和希尔伯特变换器的过渡带都非常狭窄，通带波纹均能达到指标要求，从图上看是一条平滑的直线，实际的阻带容度更略微高于所要求的容度，而且呈现出等波纹的特性（使用窗函数法设计出的波形在阻带呈现下滑趋势）。在尽量使函数波形接近于理想波形的同时，我也兼顾到了滤波器的阶数，实践可以证明，用等波纹最佳逼近法设计出的滤波器的阶数明显低于窗口法和频率取样法。

参考文献

[1] 楼顺天，李博菡.基于MATLAB的系统分析与设计[M].西安：西安电子科技大学出版社，2001，9. 200-203.

[2] 维纳.K.恩格尔，约翰.G.普罗克斯.数字信号处理—使用MATLAB[M].西安：西安交通大学出版社，2002，6. 264-280.

[3] 楼顺天，陈生潭，雷虎民.MATLAB 5.X程序设计语言[M].西安：西安电子科技大学出版社，2001，3. 1-7.

[4] 程佩青.数字信号处理教程[M].第2版.北京：清华大学出版社，2003,11. 202-226.368-380.

[5] 王超，沈伯弘.一种高速FIR数字滤波器的实现方法[J].无线电工程，1998.28 （6）.

[6] 田新广，张尔扬.一种基于窗函数法设计FIR数字滤波器的新算法[J]. 无线电工程，1999，32（8）.

[7] 成峰，周宝琨，程利青等.基于MATLAB的1/3倍频程FIR数字滤波器的设计[J].福州大学学报，2003，4.31(2).

[8] 黎雄，张学智. FIR数字滤波器的最优化设计及其MATLAB实现[J].信息技术，2004，10. 28 (10).

[9] K. Steiglitz, T. W. Parks, and J. F. Kaiser, METEOR: A Constraint-based FIR filter design

program, IEEE Trans. Signal Process., vol. 40(8), 1901–1909, 1992.

[10] A. W. Potchinkov, Design of optimal linear phase FIR filters by a semi-infinite programming

technique, Signal Process., vol. 58(2), 165–180, 1997.

[11] X. P. Lai, Design of a class of FIR filters with frequency equation constraints, Circuits Systems

Signal Process., vol. 21(2), 181–193, 2002.

[12] 赖晓平.ＦＩＲ数字滤波器的递推最小二乘设计算法[Ｊ]. 信号处理,1999,15(3).

[13] PARKS TW, MCCLELLAN J H. Chebyshev approximation for nonrecursive digital filters with linear phase [J ] . IEEE Trans Circuit Theory , 1972 , 19(2) : 1892194.

[14] 郭孔辉,轧浩.四轮转向的控制方法的发展.中国机械工程,1998，9（5）. 

[15] 闫晓艳,傅丰林,陈 健等. FIR数字滤波器的设计及其在MATLAB中的仿真实现[J].2004.25(6)

[16] 梁　静. 数字ＦＩＲ滤波器的ＭＡＴＬＡＢ设计和仿真[J].2003,13(6).

致         谢

值此论文完成之际，感谢学校图书馆和电子阅览室提供的珍贵资料；感谢信息工程学院全体老师，尤其是我的导师周小薇老师，在设计过程中得到了她的大力支持，给我提出了许多宝贵的建议和资料,同时也要感谢我的同学，没有他们的帮助，我同样无法完成本次设计。当然，由于时间仓促，本人水平所限，论文中肯定还有一些不足之处，恳请老师和同学们加以指正。

附录

1．符号说明

本程序采用国际上惯用的像形符号，例如用w或W代替；，，，在程序中分别用Wp,Ws,Rp,As代替。

在每个程序中都将直接调用一个名为freqz_m的函数，MATLAB本身并没有此函数，故首先须对其进行定义，定义的程序如下：

function [db,mag,pha,grd,w]=freqz_m(b,a);

%Modified version of freqz subroutine 

%[db,mag,pha,grd,w]=freqz_m(b,a);

%mag=absolute magnitude computed over 0 to pi radians

%pha=Phase response in radians over 0 to pi radians

%grd=Group delay over 0 to pi radians

%w=501 frequency samples between 0 to pi radians

%b=numeratou polynomial of H(z)  (for FIR:b=h)

%a=demoninatou polynomial of H(z) (for FIR:a=[1])

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501))';w=(w(1:1:501))';

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

2．设计原程序

低通滤波器

wp=0.2*pi; ws=0.3*pi; Rp=0.25; As=50;

delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1);

delta2=(1+delta1)*(10^(-As/20));

deltaH=max(delta1,delta2); deltaL=min(delta1,delta2);

weights=[delta2/delta1 1];

deltaf=(ws-wp)/(2*pi);

M=ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1);

M=M+4;

f=[0 wp/pi ws/pi 1];

m=[1 1 0 0];

h=remez(M-1,f,m,weights);

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

wsi=ws/delta_w+1;wpi=wp/delta_w;

Asd=-max(db(wsi:1:501))

M

%

%Plots

%

n=[0:1:M-1];

subplot(2,2,1),stem(n,h),axis([0,50,-0.1,0.3]),xlabel('n'),ylabel('h(n)'),title('Actual Impulse Response'),

set(gca,'XTickMode','manual','Xtick',[0,M]),set(gca,'YTickMode','manual','YTick',[-0.1:0.1:0.3]),hold;

subplot(2,2,2),plot(w/pi,db),axis([0,1,-80,10]),xlabel('\omega(单位为\pi)'),ylabel('dB'),title('Magnitude Response in dB'),

set(gca,'XTickMode','manual','Xtick',[0,0.2,0.3,1]),set(gca,'YTickMode','manual','YTick',[-50,0])

subplot(2,2,3),plot(w/pi,mag),axis([0,1,-0.1,1.1]),xlabel('\omega(单位为\pi)'),ylabel('Hr(\omega)'),title('Amplitude Response'),

set(gca,'XTickMode','manual','Xtick',[0,0.2,0.3,1]),set(gca,'YTickMode','manual','YTick',[0,1])

运行结果滤波器长度为47，阻带容度51.0857分贝。

高通滤波器

ws=0.6*pi; wp=0.75*pi; Rp=0.5; As=50;

delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1);

delta2=(1+delta1)*(10^(-As/20));

deltaH=max(delta1,delta2); deltaL=min(delta1,delta2);

weights=[1 delta2/delta1];

deltaf=(wp-ws)/(2*pi);

M=ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1);

M=2*floor(M/2)+1

M=M+2;

f=[0 ws/pi wp/pi 1];

m=[0 0 1 1];

h=remez(M-1,f,m,weights);

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

wsi=ws/delta_w;wpi=wp/delta_w;

Asd=-max(db(1:1:wsi))

M

%

%plots

n=[0:1:M-1];

subplot(2,2,1),stem(n,h),axis([0,28,-0.4,0.4]),xlabel('n'),ylabel('h(n)'),title('Actual Impulse Response'),

set(gca,'XTickMode','manual','Xtick',[0,M]),set(gca,'YTickMode','manual','YTick',[-0.4,-0.2,0,0.2,0.4]),hold

subplot(2,2,2),plot(w/pi,db),axis([0,1,-80,10]),xlabel('\omega(单位为\pi)'),ylabel('dB'),title('Magnitude Response in dB'),

set(gca,'XTickMode','manual','Xtick',[0,0.6,0.75,1]),set(gca,'YTickMode','manual','YTick',[-50,0])

subplot(2,2,3),plot(w/pi,mag),axis([0,1,-0.1,1.1]),xlabel('\omega(单位为\pi)'),ylabel('Hr(\omega)'),title('Amplitude Response'),

set(gca,'XTickMode','manual','Xtick',[0,0.6,0.75,1]),set(gca,'YTickMode','manual','YTick',[0,1])

运行结果滤波器长度为29，阻带容度50.2253分贝。

带通滤波器

ws1=0.2*pi; wp1=0.35*pi; wp2=0.65*pi; ws2=0.8*pi;

Rp=1.0; As=60;

delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1);

delta2=(1+delta1)*(10^(-As/20));

deltaH=max(delta1,delta2); deltaL=min(delta1,delta2);

weights=[1 delta2/delta1 1];

delta_f=min((ws2-wp2)/(2*pi), (wp1-ws1)/(2*pi));

M=ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*delta_f)+1);

M=M+1

f=[0 ws1/pi wp1/pi wp2/pi ws2/pi 1];

m=[0 0 1 1 0 0];

h=remez(M-1,f,m,weights);

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;

ws1i=floor(ws1/delta_w)+1; wp1i=floor(wp1/delta_w)+1;

ws2i=floor(ws2/delta_w)+1; wp2i=floor(wp2/delta_w)+1;

Asd=-max(db(1:1:ws1i))

%

% Plots

%

n=0:M-1;

subplot(2,2,1),stem(n,h),xlabel('n'),ylabel('h(n)'),title('Actual Impulse Response'),

set(gca,'XTickMode','manual','Xtick',[0,M]),set(gca,'YTickMode','manual','YTick',[-0.4:0.2:0.4]),hold;

subplot(2,2,2),plot(w/pi,db),axis([0,1,-90,10]),xlabel('\omega(单位为\pi)'),ylabel('dB'),title('Magnitude Response in dB'),

set(gca,'XTickMode','manual','Xtick',[0,0.2,0.35,0.65,0.8,1]),set(gca,'YTickMode','manual','YTick',[-60,0])

subplot(2,2,3),plot(w/pi,mag),axis([0,1,-0.1,1.1]),xlabel('\omega(单位为\pi)'),ylabel('Hr(\omega)'),title('Amplitude Response'),

set(gca,'XTickMode','manual','Xtick',[0,0.2,0.35,0.65,0.8,1]),set(gca,'YTickMode','manual','YTick',[0,1])

运行结果滤波器长度29，阻带容度61.2843分贝。

希尔伯特变换器

f=[0.05,0.95]; m=[1 1];

h=remez(50,f,m,'hilbert');

[db,mag,pha,grd,w]=freqz_m(h,[1]);

subplot(2,1,1); stem([0:50],h); title('Impulse Response');

axis([0,50,-0.8,0.8]);xlabel('n');ylabel('h(n)');

set(gca,'XTickMode','manual','XTick',[0,50])

set(gca,'YTickMode','manual','YTick',[-0.8:0.2:0.8]);

subplot(2,1,2); plot(w/pi,mag); title('magnitude Response')

xlabel('frequency in pi units');ylabel('H')

set(gca,'XTickMode','manual','XTick',[0,f,1])

set(gca,'YTickMode','manual','YTick',[0,1]); grid

运行结果长度51，处于频带之内。

首页 上一页 1 2 3 下一页 尾页 3/3/3

相关论文