信号处理工程应用训练(指导书).docx

上传人:小飞机 文档编号:1851847 上传时间:2022-12-21 格式:DOCX 页数:52 大小:295.13KB
返回 下载 相关 举报
信号处理工程应用训练(指导书).docx_第1页
第1页 / 共52页
信号处理工程应用训练(指导书).docx_第2页
第2页 / 共52页
信号处理工程应用训练(指导书).docx_第3页
第3页 / 共52页
信号处理工程应用训练(指导书).docx_第4页
第4页 / 共52页
信号处理工程应用训练(指导书).docx_第5页
第5页 / 共52页
点击查看更多>>
资源描述

《信号处理工程应用训练(指导书).docx》由会员分享,可在线阅读,更多相关《信号处理工程应用训练(指导书).docx(52页珍藏版)》请在三一办公上搜索。

1、信号处理工程训练课训练一 信号与系统函数编程训练目的 1学会将信号与系统函数转变成计算机程序。 2基本掌握将数学函数转变为程序函数的技巧与规范。 3了解理论函数与程序函数的差异。初步认识计算机适用范围。 训练介绍 1数学函数转化问题把根据数学函数编写的C函数子程序称为程序函数。数学函数与程序不可能完全一致。一是计算机运算都有一个范围,所做运算超出范围便会出错;二是因为计算机不能做除零运算,这会产生一除法错,理论函数无此限制。所以要求在编写程序函数时一定要结合实际应用情形来确定如何编写,不能简单照搬数学函数。三是程序函数不象数学函数那样易于进行代数运算或者具有某种运算性质,例如理论上的冲击函数,

2、则不易编写对应的函数子程序,所以数学函数并不能全由计算机的程序函数完全实现。一般在将一数学函数转变为一计算机上程序函数时,要具体情况具体处理。编写程序函数有一些规范和注意事项:(1)数学函数当中若有除法运算,需仔细函数奇异值的处理,须通过程序中的判断和特殊处理使程序函数返回正确值。(2) 数学函数中跳变点的极限值,常取左右极限的均值,程序函数中以右极限作为函数的取值。若特殊需要,须与数学函数完全一致,则仍按数学函数规定取值。(3) 所有函数子程序的输入与输出参量尽量规定为double型,建议不用float型,这是出于规范考虑。(4) 所有程序函数的输入输出参量声明时写成如下形式: Double

3、 function(Type out1,Type out2,. Type in1,Type in2,.) Double function(Type out1,Type out2,. Type io1,Type io2,. Type in1,Type in2,.)即,输出变量占一行,输入输出变量占一行,输入变量占一行。输入变量的第一个参量为主变量。(5) 尽量减少函数变量个数,例如sin(wt)有两个参数,编程只需实现sin(x)。(6) 每个函数子程序须有适当文字注释,注释的内容包括索引号,对应的理论函数,编者姓名及日期,函数的功能定义域值域,使用举例等。说明应简洁清楚,以备能长期正确使用。(

4、7) 程序函数块内的小块以一空行进行分割,程序函数体之间,以2、3空行行分割。组织一个函数库文件时应将功能,特征相近的函数子程序归在一起。各分类块间应有适当的注释说明。 2以下以单位阶跃U(t)、方波和函数三种信号函数为例进行编程示范:训练内容 0 1、斜变函数R(t)= t,t0 2、锯齿波:f(t)=t / T,0t。训练步骤 1依训练内容在一个文件中编写好三个函数子程序,并依要求进行注释,做成一个库文件。 2编写一个有main( )函数的可执行文件,在此执行文件中调用库文件里的函数子程序进行计算,计算结果由计算机屏幕输出,结合手算验证,多取不同情况的特殊值,保证程序函数正确。问题讨论 1

5、举例说明不宜用程序函数表达的理论函数,并总结这些函数的特点。2讨论理论方式、计算机方式在处理实际物理的信号与系统问题时的异同。3对规范化程序方法的初步认识。训练二 图形显示与观察训练目的 1掌握基本的计算机作图方法。 2掌握常用的信号与系统图形观察。 3熟悉一些显示器编程的基本概念。 训练介绍 1、计算机绘图基础对于计算机屏幕,机器本身按物理坐标绘图显示,屏幕左上角的物理坐标是(0,0),右下角的坐标是(横向分辨率-1,纵向分辨率-1)。Windows操作系统是在屏幕上开若干窗口,有客户区的窗口用户可以在客户区作图,窗口客户区的左上角的相对物理坐标是(0,0),右下角的相对物理坐标是(w-1,

6、h-1),w是宽,h是高,用户工程使用的是自己定义的用户坐标,用户坐标到窗口的相对物理坐标有一个转换。如图2-1所示: (图2-1)用户坐标可以是三维坐标,同样是转换到2维屏幕窗口坐标。工程中常用不同的窗口模拟一些仪器,如用示波器扫频仪频谱仪矢量分析仪等等,也可用不同的坐标系统和作图函数工具研究分析工程中的信号和系统。对窗口的作图由不同的图形驱动函数实现。基本的作图函数有SetPixel()、MoveToEx()、LineTo(),其它作图函数可由基本函数构造出来。工程训练中,对2维用户坐标窗口作图函数,函数名后有一个2,如linto2(),3维的作图函数后面有一个3。 不同窗口不同坐标体系的

7、作图函数众多,主要的是掌握基本,理解其它或构建其它。训练内容 1、用plotxy2( )绘制余弦,正切,和e为底的指数函数; 2用plotxy2( )绘制训练一中的斜变函数和锯齿波函数,用plotxyz3()绘制sinC()函数。 3、练习物理坐标下的作图,绘制圆,三角形和扇型图(参数自定)5北京工业大学电控学院信号处理工程训练课训练三 波形合成训练目的 1学会用计算机合成波形。 2学习将理论知识同用计算机相结合的方法。 训练介绍波形合成有许多运用,如电子琴、信号源、计算机里声卡,合成法使得信号、声音、图象的产生更加容易,使用更加方便,计算机研究信号波形合成,可直观看出合成的具体效果,这对我们

8、进行科学研究和设计新仪器新设备都有极大帮助。训练理论,核心是一种函数逼近,通过合成信号g(t)希望能与理想信号f(t)差别很小,以使实际使用当中可用合成信号代替理想信号。用数学形式表示为:g(t)=c1g1(t)+c2g2(t)+gn(t),(g1(t)g2(t) . gn(t)为一函数集)在(3-1)条件下的c1,可得到f(t)与g(t)的最佳匹配。由此可找到函数f(t)在不同坐标函数集上的分解合成形式。 其中由正交正弦函数集表示的傅立叶级数为: 式中 有关系函数介绍:积分工具函数为宏函数_Trapz(value,a,b,N,fx,x)将自变量为x的f(x)在a,b区间划分成进行梯形积分,返

9、回的积分值,注意x为形式上参数,用于指明fx函数表达式中的自变量,无须用前声明。double FouSer(double t,double *T, double a,double b,int *N);用于求(3-2)式的值。t:为(3-2)中的自变量;T:为信号与系统函数的周期;a,b:为(3-2)中系数构成的数组,数组长度为N。 一个周期为1的函数f(t)的正交沃尔什级数表示为: 其中: 这里kr是k的二进制各位数字,p是二进制位数。 Walsh函数只取+1,-1,例如Wal(1,t)=sgn(cospt)在0=t1/2时为+1,在1/2=t 4key add dx=(b-a)/n; _ ,

10、 ”四键为增,“-,”四键为减),x0变化区间为a,b,变化步长为(b-a)/n。按ESC键函数返回0值,其他键返回1。int getkey(void);预定义在xxgc.h;无返屏地返回一次普通按键或附加按键的值。示例对一周期信号扩展到整个t轴。分别产生上述三种失真,编程进行模拟(注意先确定好不失真的输入信号是什么)。训练内容 1以训练三的k级合成三角波和方波为不失真的原始信号(k由自选),总结以下情况失真的特点。(a) 动态观察波形在某一谐波处的单纯幅度、相位、频率失真。(b) 研究某种线性组合的幅度、相位、频率失真。(c) 研究同时发生有幅度相位、相位频率、幅度频率、幅度相位频率的失真。

11、问题讨论 1总结各种失真特点,讨论失真可由示波器判别的方法。训练五 卷积运算训练目的 1学习用数值计算方法实现卷积积分。 2学习构造通用的函数工具。训练介绍任何实际物理系统与信号都可认为有起始点,所以普通数学意义的卷积可以是:用编程实现卷积运算,简单地,可直接套用(5-2)式的矩形积分运算式,或使用积分工具宏函数_Trapz()完成积分。训练内容 1用方法一二实现下面的卷积分,并作出卷积图形:训练六 卷积图解训练目的 1学习制作简单演示软件。 2通过动态图形观察进一步理解卷积概念。 训练介绍 卷积是种形象叫法,这是因为两函数作卷积分,好象这两函数互相卷绕在一起一样,要展示这一个过程,可用作图方

12、法将其形象表现出来。让卷积的两函数之一反转并从某一时刻开始平移,两函数图形将会慢慢重迭,重迭部分面积就是运动某一时刻的卷积值。实际上,从图形动态移动可以清楚看出,两图形并为发生卷绕,只是对褶后的平移,所以有的书称褶积。卷积图解过程是先将一函数图形沿纵轴反转,然后从左到右平移,反转图形在不同平移位置与另一函数图形重叠部分的乘积代表不同平移时刻卷积值。在主函数中要实现动态卷积演示效果,可用控件函数instKeyCtr()与keyCtr()控制反转函数平移位置。反转平移的函数卷积结果未反转的函数 (图6-1)训练内容 自选两个f(t),h(t)作动态卷积分观察。训练七 周期信号频谱观察训练目的 1学

13、会周期信号频谱绘制方法。 2学会数值积分函数的使用及技巧。 训练介绍 求解周期信号傅立叶级数的公式为:T为信号周期,w=2p/T 。单边频谱的振幅为:A0=a0, n=123. 相位谱为:Qn=arctgbn/an , n=123. (7-2)写成复数级数形式为: (7-3)根据求出的an 、bn,构造出复数级数系数,以n为自变量,在3维坐标中用绘图工具函数plotgri()画出所有复数系数,示例以窗函数为例。训练内容 1动态观察三角波和升半正弦波的幅度相位频谱特性。训练八 付氏变换研究训练目的 1学会用计算机绘制函数付氏变换后的频谱图形。 2学习特殊函数积分的转化方式。 训练介绍付氏变换在理

14、论上由下面变换对实现:这里仅考虑(8-1)式,被积函数为一复数,为了能利用实数积分函数子程序iTrapz(),先将(8-1)式作如下分解:计算机一般不能处理(8-3)中无穷区间的积分,只能处理被积函数在有限区间取值有限的积分,对于无穷区间积分,假设函数自变量在趋于无穷当中,函数值也趋于0,可先取有限区间对函数积分,然后将积分区间扩大一倍,再由计算机算出积分结果,比较前后两次结果,继续扩大积分区间,直到两次积分结果相差很小可以忽略,或者没有差别,那么,有理由认为,再增加积分区间,也不会对积分值有多大贡献,所以这个足够大区间上的积分值,就可以用来表示无穷区间上的积分。考虑函数的傅氏变换,理论计算得

15、根据(8-4)编写程序函数fetut()作为标准参考,以对照用计算机方法绘制的f(t)的傅立叶变换频谱曲线正确与否。 训练内容 1、用计算机作出以下函数的傅氏频谱: (a) 半余弦脉冲 (b) 三角脉冲 (c) 高斯脉冲 (d) 抽样函数 2、构造傅氏反变换积分工具。3、作上面函数在二维用户坐标下的幅度相位频谱。19北京工业大学电控学院信号处理工程训练课训练九 系统的幅频特性显示与观察训练目的 学会用计算机作系统函数的幅频相频特性曲线。 训练介绍在正弦稳态下,系统的响应向量与激励向量之比称为系统函数。其一般表达式为: N(s) a0sn + a1sn-1 + . + anH(s) = = (9

16、-1) D(s) b0sm + b1sm-1 + . + bm其中:s = jw 。系统函数的模|H(jw)|与频率关系称为“幅频特性”。系统函数的幅角f(w)表示了响应与激励的相位差,其随频率的关系称为“相频特性”。 示例以为例,编程绘制其幅频、相频特性。 训练内容 1用计算机作出以下系统传递函数的幅频与相频特性曲线: 10s2(1) H(s) = s2+2s+10 200(2) H(s) = (s+1000)(s-2000) 1+s/10(3) H(s) = 10 (1+s)(1+s/100)2作系统的波特图(横坐标取10为底的对数,纵坐标20log)。3、 用plotwri()绘幅度相位

17、图。训练十 状态方程的数值计算训练目的 掌握由状态变量法求解系统响应的计算方法。 训练介绍当采用状态变量描述系统时,系统的状态方程与输出方程可用下面矩阵形式表示: (10-1) (10-2)求解(10-1)的一阶微分方程,常用的数值方法有:尤拉法梯形法和龙格库塔法。此次训练采用最简单的尤拉法,考虑下面方程: (10-3)具有初始条件:x(t0)=x0从初始值x(t0)开始,按步计算x(t1)x(t2).x(tn),用第一差分代替导数,得或者写成为方便起见,对于所有k的递增常取相同的量,用h表示,称为步长。令tk=t0+hk及x(tk)=xk,则有xk+1=xk+f(xk,tk)h (10-4)

18、此即为尤拉公式。 尤拉近似的要点,在于把每一时间间隔之内的函数的导数看作常数,且等于间隔起始时的导数值。 为说明这种算法,现设有二阶线性时不变系统的两个状态微分方程和一个输出方程如下: y=c11x1+c12x2+de尤拉近似法的计算步骤如下: 训练内容 1、设已知系统的状态方程和输出方程为: 1试用尤拉近似法求解,要求输出x1x2和y在02秒内的值及波形。2、 编写通用的尤拉法求解的函数工具。训练十一 DFT性质研究训练目的 1了解DFT正反变换的计算关系。 2了解DFT实虚部的对应关系。 3了解x(n)补零后,对DFT后X(k)产生的影响,或X(k)补零后对x(n)的影响。 训练内容 1自

19、编一DFT计算子程序,要求该子程序具正反变换的双重功能,计算用复数实现。 形式为:void DFT(COMPLEX outArray, COMPLEX inArray,int n,BOOL flag) 2令x(n)=n+1,n=0,1,.,7,显然,x(n)是x(t)=t+1得到的,是一斜坡信号,由它试验所编程序是否正确(即:对x(n)做DFT得X(k),再由X(k)做反变换,应得同样的x(n)。验证能量守恒关系(即Parseval定理): 3将x(n)补n个零,即求X(k)和分析X(k)的关系,并给出能量关系。 4对x(n)=n+1,n=0,1,.,7,作插值,按 x(2n)=x(n),n=

20、0,1,.,7 x(2n+1)=x(n)的某种线性组合,n=0,1,.,7试分析应为何实现,并实现之。 5令x(2n)=x(n)=n+1,n=0,1,.,7x(n)=0,n=0,1,.,7x(n)实为x(n)在相邻点间补零而成。求X(k)。注:对训练345请在报告中给出理论推导。训练十二 DFT及抽样定理研究训练目的 了解线谱的特性DFT的性质及对正弦信号的抽样方法。训练介绍抽样定理指出,若信号x(t)的最高频率为fc,当抽样频率fs2fc时,可由抽样信号x(nTs)完全恢复原信号x(t)。该定理对正弦信号有其特殊性。我们知道,一个正弦信号x(t)=sin(2pfct),其频谱在fc处,是一d

21、函数,这样的谱称之为线谱。具线谱特性的典型信号是正余弦函数,它们广泛应用于信息技术领域,因此研究正弦信号的抽样是正确进行其它信号处理工作的基础。训练内容 1给定信号 x(t)=sin(2pfct),fc=50,N=264(1) 按下面给的fs对x(t)抽样,求时域能量(2) 对x(n)做DFT,求出X(k),并求出频域能量 记 在50hz处的频谱为X50,若抽样频率正确,即无泄漏,则 Et=Ef=2|X50|2/N (注:这里用到多项式插值函数inspoly()(3) 观察有无泄漏,若有,是多少?并分析产生泄漏的原因。(4) 对X(k)做反变换,再求x(n),观察恢复后的x(n)与原x(n)有

22、无区别请按下面抽样频率,重复上面4个步骤:(1) fs=100hz(2) fs=110hz(3) fs=200hz(4) fs=230hz(5) fs=250hz 2仍令fs=250hz(1) 取N=250点,用DFT直接求X(k);k=0,1,.,249(2) 取N=250点,补零至N=256,用FFT求X(k)按训练1的四项要求进行,通过训练分析,在什么条件,可保证对正弦信号的抽样无泄漏。 3给定正弦信号的频率fc,FFT的长度N=2m,若要不泄漏,抽样频率fs应如何选择52北京工业大学电控学院数字信号处理训练十三 数字滤波器制作训练目的 1学习制作简单的数字滤波器。训练介绍由于数字滤波器

23、有着体积小运用灵活精度高可靠性好等优点,已在许多领域获得应用,如今被广泛用于数字电视,数字通讯,计算机控制,信息分析等诸多方面,成为数字革命的象征。数字滤波器不同与传统的LC陶瓷模拟等滤波器,是由于它是在人们深入分析了传统滤波器本质之后,实现的一种同功同效。其结构如下: 图( 13-1)它将A/D得到数字信号经过数学运算,运算结果由D/A输出。 数字滤波器系统函数的形式一般为: (13-1) 将其写成差分方程: ( 13-2) 若N=1则称为一阶数字滤波器: y(n)=b0x(n)+b1x(n-1)+a1y(n-1) ( 13-3)它的一种结构实现如下图所示: 图(11.2)一阶数字滤波器结构

24、图 若N=2则称二阶数字滤波器: ( 13-4)它的构成方式有多种,一种称为直接型结构,如图(11.3)所示。 图( 13-3)注:在结构图中,以移位算子z-1代表延时,以箭头旁的加权系数代表乘以常系数的乘法运算,若系数为1,则可不标,若某项系数为0,则箭头应断开。 下面是图( 13-2)的实现流程图(P为序列中数据的个数)。 图(11.4) 一阶数字滤波器流程图 训练内容 1设一阶数字滤波器的系统函数为: 系统输入为以下三种数字序列:(1) 矩形序列 图( 13-5)(2)三角序列 图( 13-6)(3)正弦序列 图( 13-7)绘出三种数字序列信号经一阶数字滤波器后的输出序列。 2设二阶数

25、字滤波器的系统函数为: 其输入为单位阶跃信号 图( 13-8) 3、编写通用数字滤波器函数工具。训练二十 IIR数字滤波器设计与实现训练目的掌握用双线性变换法设计IIR数字滤波器的方法;掌握由低通到高通带通带阻的频率变换方法;掌握Butterworth滤波器设计的原理。训练介绍1、标准形式ButterworthAF到DF的转换a) 把Butterworth 低通AF在db1分贝衰减处的频率规定为1,频率归一的Butterworth LP-AF形式为: ( 14-1)则, (14-2) 简单标准形式为: (14-3)简单标准形式转变成普通频率归一AF的关系为:b) 由频率归一的Butterwor

26、th低通AF转换成非归一的其他类型滤波器 图(14-1) 频率归一的butterworth LP AF幅频特性LP到LP: ,1W1 ( 14-4) 图(14-2)LP到HP: ,-1W2 ( 14-5) 图(14-3)LP到BP:,1W2, -1W3 ( 14-6) 图(14-4)LP到BS:,1W1, -1W4 ( 14-7) 图(14-5)c) 由AF到DF的变换,即双线性变换为: ( 14-8)根据双线性变换公式: 知DF频率f与对应AF角频率W的关系为: , (14-9) 将a)、b)、c)三步合成一步, 可得简单标准形式的Butterworth LP AF到DF的代换关系:DF为L

27、P: (14-10)DF为HP: ( 14-11)DF为BP:, ( 14-12)DF为BS:, ( 14-13)在C语言编程时,把( 14-10)( 14-11)( 14-12)( 14-13)统一看作 ( 14-14)形式的有理分式,用C语言函数btwC23( )返回分式分子分母的系数,double btwC23(double23,int,int,double,double,double,double,.);double btwC23(c,bandType,N,db1,fs,f1,f2,f3,f4)int bandType,N;double c23,db1,fs,f1,f2,f3,f4;其

28、中输出参数:c数组,在函数运行后装填的是( 14-14)式的系数;函数的返回值,是对应DF的AF的归一用频率。输入参数:bandType为通带类型,可选用枚举值LOWPASS、HIGHPASS、BANDPASS、BANDSTOP; N为LP Butterworth的阶次; fs为DF的抽样频率,fs=1/T; db1,f1,f2,f3,f4规定为(图 14-6)示 20lg|H(z)| 20lg|H(z)| 20lg|H(z)| 20lg|H(z)| f1 f2 f f1 f2 f f1 f2 f3 f4 f f1 f2 f3 f4 f-db1 -db1 -db1 -db1-db2 -db2

29、-db2 -db2LOWPASS HIGHPASS BANDPASS BANDSTOP(图 14-6)函数参数调用示意()2、计算Butterworth LP AF的阶次N:根据( 14-2)式,不同类型DF对应的Butterworth LP AF都有 ( 14-15)由LP DF参数求N,阻带内: ( 14-16)联立( 14-15)、( 14-16)解得: (14-17)由HP DF参数求N,阻带内: (14-18)联立( 14-15)、( 14-18)解得: ( 14-19)( 14-17)和( 14-19)表明LP DF和HP DF求阶次N的形式一样。由BP DF参数求N,设频率映射,

30、阻带内: ( 14-20) ( 14-21)其中由( 14-20)( 14-21)和( 14-15)联立解得: 和 ( 14-22)选择和小的N作为结果。由BS DF参数求N,设频率映射,阻带内: ( 14-23) ( 14-24)其中由( 14-23)( 14-24)和( 14-15)联立解得: 和 ( 14-25)选择和小的N作为结果。将( 14-17)( 14-19)( 14-22)( 14-25)的计算过程编制成C语言函数,阶次N可由C函数btwOrder( )的返回值得到,函数的声明为:int btwOrder(int,double,double,double,double,doub

31、le,)int btwOrder(bandType,db1,db2,fs,f1,f2,f3,f4)int bandType;double db1,db2,fs,f1,f2,f3,f4;其中输入参数: bandType为滤波器的类型,可选用枚举值LOWPASS,HIGHPASS,BANDPASS,BANDSTOP;db1为通带衰减分贝值;db2为阻带衰减分贝值;fs为数字滤波器的抽样频率;f1,f2,f3,f4如(图 14-6)所示。3根据阶次N设计简单标准形式的Butterworth LP AF简单标准形式的Butterworth LP AF为 ( 14-26)令分母多项式为0,求其根得: (

32、 14-27)所有左半平面的极点为满足的点,即,也即0kN时的极点为稳定极点,且Sk与SN-k-1为一对共轭极点,所以,简单标准形式的二阶Butterworth LP AF的转移函数可写为: ( 14-28) 在N为偶数, ( 14-29a) N为奇数时, ( 14-29b)在C语言编程时,把( 14-29a)( 14-29b)统一看作 ,L=N+1/2取整 ( 14-30)的形式,用C语言的二维数组b表示butterworth LP AF 系统函数的分母多项式系数,C函数声明为:int btwAf(double b2,int N)函数的输入参数:N为butterworth LP的阶次;输出参

33、数:b为预先声明的一个长度为L的二维数组,函数返回时b数组内填好( 14-30)的系数; 函数的返回值为L。 4、简单标准形式butterworth LP AF向DF转换的程序自动实现把( 14-14)代入( 14-30)得: ( 14-31)为书写方便,用hijk表示C语言形式hijk, , , , , 实现此一过程的函数:void btwAf2Df(double h25,int L, double b2,double c23)5、DF幅频特性的绘制将z=ejT=ej2fT代入(14-31),则函数|H(z)|变成以f为自变量的函数,能根据系统函数数组h,抽样频率fs,频率f返回20lg|H

34、(ej2fT)|的值的C函数声明为:double btw20lgHz(double f,double fs,double h25,int L);以20lg|H(ej2fT)|为纵轴,以f为横轴,用ploxy2( )绘制曲线,图形输出结果可参照图(14-6)。6、用DF系统函数对信号进行滤波训练内容 1(1)令:f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。 设计LP DF(2)令:f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。 设计HP DF(1) 令:f1=200hz,f2=300hz,f3=400hz,f4=500hz,fs=2000hz,a1=3dB, a2=40dB。设计BF DF(2) 参数同3,设计BS DF注:每种滤波器都应给出(1)转移函数(2)幅频相频曲线2画出整个设计过程的程序流程图,说明各部分功能。3、实现滤

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 生活休闲 > 在线阅读


备案号:宁ICP备20000045号-2

经营许可证:宁B2-20210002

宁公网安备 64010402000987号