IIR带通滤波器语音去噪要点.pdf

上传人:tbuqq 文档编号:5197176 上传时间:2020-02-19 格式:PDF 页数:20 大小:306.14KB
返回 下载 相关 举报
IIR带通滤波器语音去噪要点.pdf_第1页
第1页 / 共20页
IIR带通滤波器语音去噪要点.pdf_第2页
第2页 / 共20页
IIR带通滤波器语音去噪要点.pdf_第3页
第3页 / 共20页
IIR带通滤波器语音去噪要点.pdf_第4页
第4页 / 共20页
IIR带通滤波器语音去噪要点.pdf_第5页
第5页 / 共20页
点击查看更多>>
资源描述

《IIR带通滤波器语音去噪要点.pdf》由会员分享,可在线阅读,更多相关《IIR带通滤波器语音去噪要点.pdf(20页珍藏版)》请在三一文库上搜索。

1、摘要 语音信号滤波处理时研究用数字信号处理技术和语音学知识对语音信号进 行处理的新兴学科,是目前发展最为迅速的信息科学研究领域的核心技术之一, 通过语音传递信息是人类交流信息最自然、最有效、最方便的手段。 本次主要通 过录制一段语音,对其进行时域、频谱分析,并利用matlab 的信号处理工具箱 对语音进行加噪然后再用IIR 数字带通滤波器滤除噪声, 最后对比滤波前后的语 音信号的时域、频域特性。 关键字: IIR ;双线性变换;模拟低通滤波器;切比雪夫;MATLAB 目录 目录 2 前言 1 一设计原理 . 2 1. 数字滤波器简介. 2 2.IIR数字滤波器的设计原理 3 3.IIR滤波器的

2、特点 3 二 IIR 数字滤波器的设计方法 4 1. 模拟滤波器 . 4 2. 双线性变换法. 6 三 IIR 数字滤波器设计过程 9 1. 设计步骤 . 9 2. 音频信号部分程序. 10 3. 程序流程图 . 10 . 仿真结果 . 11 总结 13 致谢 14 参考文献 . 15 附录: 16 1 前言 通过语音传递信息室人类最重要、最有效、最常用和最方便的交换信息的形 式。语音是人类特有的功能, 声音是人类最常用的工具, 是相传递信息的最重要 的手段。因此,语音信号是人类构成思想疏通和感情交流的最重要的途径之一。 并且,由于语言和语音与人的智力活动密切相关,与社会文化和进步紧密相连,

3、所以它具有最大的信息容量和最高的智能水平。现在,人类已开始进入了信息化 时代,用现代手段研究语音信号,使人们能够更加有效地产生、传输、存储、获 取和应用信息, 这对于促进社会的发展具有十分重要的意义。让计算机能听懂人 类的语言,是人类自计算机诞生以来梦寐以求的想法。 随着计算机越来越向便携化方向发展,随着计算环境的日趋复杂化, 人类越 来越迫切要求拜托键盘的束缚而带至以语音输入这样便于使用的、自然的、人性 化的输入方法。 作为高科技应用领域的热点, 语音信号的采集和分析从理论的研 究到产品的开发已走过了几十个春秋并且取得了长足的进步。它正在直接与办 公、交通、金融、公安、商业、旅游等行业的语音

4、咨询与管理,工业生产部门的 语音声控, 电话、电信系统的自动拨号、 辅助控制和查询以及医疗卫生和福利事 业的生活支援系统等各种实际应用领域相接轨,并且有望成为下一代操作系统的 和应用程序的用户界面。 可见,语音信号采集与分析的研究将是一项极具市场价 值和挑战性的工作。 我们今天进行这一领域的研究与开发就是要让语音信号处理 技术走向人们的日常生活当中,并不断朝更高目标而努力。 语音信号采集与分析之所以能够那样长期地、深深地吸引广大科学工作者对 其进行研究和探讨, 除了它的实用性之外, 另一个重要原因是, 它始终与当时信 息科学中最活跃的前沿的学科保持密切连续,并且一起发展。 语音信号的采集和 分

5、析是以语音语言学和数字信号处理为基础而形成的一门设计面很广的综合性 学科,与心理、生理学、计算机科学、通信与信息科学以及模式识别和人工智能 等学科都有着非常密切的关系。 对语音信号的采集和分析的研究一直是数字信号 处理技术发展的重要推动力量。 因为许多处理的新方法的提出, 首先是语音信号 处理中获得成功,然后在推广到其他领域。 2 一设计原理 IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法 来进行设计,通常采用的模拟滤波器原型有butterworth函数、chebyshev 函数、 Bessel 函数、椭圆滤波器函数等。 1. 数字滤波器简介 数字滤波器是一种用来过滤时间

6、离散信号的数字系统,通过对抽样数据进行 数学处理来达到频域滤波的目的。可以设计系统的频域响应, 让它满足一定的要 求,从而对通过该系统的信号的某些特定的频率成分进行过滤,这就是滤波器的 基本原理。 如果系统是一个连续系统, 则滤波器称为模拟滤波器。 如果系统是一 个离散系统,则滤波器称为数字滤波器。数字滤波器的技术指标如图1-1 所示: 图 1-1 IIR数字滤波器的技术指标 p: 通带截止频率 s: 阻带截止频率 s: 阻带波动p: 通带波动 通带衰减( dB) : 阻带衰减( dB): 1 0 通带 过渡带 阻带 ) ( j e H ps s p 1 )1lg(20 pp A ss lg2

7、0A 3 信号通过线性系统后, 其输出就是输入信号和系统冲击响应的卷积。除此之 外,输出信号的波形将不同于输入信号的波形。从频域分析来看, 信号通过线性 系统后,输出信号的频谱将是输入信号的频谱与系统传递函数的乘积。除非系统 函数是常数,否则,输出信号的频谱将不同于输入信号的频谱。 2.IIR数字滤波器的设计原理 (1)按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指 标; (2)根据模拟滤波器的技术指标设计为相应的模拟低通滤波器; (3)根据双线性不变法把模拟滤波器转换为数字滤波器; (4)根据要设计的带通滤波器,首先把它们的技术指标转化为模拟低通滤波器 的技术指标, 设计为相

8、应的模拟低通滤波器, 最后通过频率转化的方法来得到所 要的滤波器如下图1-2 所示: zss spp 双线性不变法设计模拟滤波器 、 频率变化 、 图 1-2 IIR数字滤波器的设计原理 3.IIR滤波器的特点 (1)IIR 数字滤波器的系统函数可以写成封闭函数的形式。 (2)IIR 数字滤波器采用递归型结构, 及结构上带有反馈环路。 IIR 滤波器运算 结构通常由延时、 乘以系数和相加的基本运算组成,可以组合成直接型、 并联型、 级联型、正准型四种结构形式,都具有反馈回路,由于运算中的舍入处理,使误 差不断累积,有时会产生微弱的几声震荡。 (3)IIR数字滤波器在设计上可以借助成熟的模拟滤波

9、器的成果,如:巴特沃 斯、切比雪夫和椭圆滤波器等, 有现成的设计数据和图表可查,其设计工作量比 较小,对计算工具的要求不高。 在设计一个 IIR 数字滤波器时, 我们根据指标先 写出模拟滤波器的公式, 然后通过一定的变化, 将模拟滤波器的公式转化为数字 滤波器的公式。 (4)IIR数字滤波器的相位特性不好控制,对相位要求较高时,需加相位校准 网络。 4 二 IIR 数字滤波器的设计方法 IIR数字滤波器是一种离散时间系统,其系统函数为公式(2-1) : N k p M r r N k k k M k k k zz zz k za zb ZH 1 1 1 0 )( )( 1 )( (2-1) (

10、1)当 M=N,H(z):N 阶 IIR 系统+(M-N )阶 FIR 系统。 (2)以上两种表示等价,部分分式形式和零极点增益形式。 (3)IIR 系统的逼近,就是找到滤波器的系数 k a , k b,或者是系统的零极点和 增益( z,p,k ) 。 1. 模拟滤波器 (一)技术指标如图2-1 所示: 图 2-1 模拟滤波器的技术指标 p:通带截止频率 s;阻带截止频率 p : 通带波动 s : 阻带波动 通带衰减( dB) : |H( j)| 1 0 通带 过渡带 阻带 ps s p 1 )1lg(20 pp A 5 阻带衰减( dB): 滤波器的衰减函数:jwwlg20dB (二)切比雪

11、夫 I 型的模拟低通滤波器的幅度响应如下图2-2 所示: 图 2-2 切比雪夫I 型的模拟低通滤波器的幅度响应 幅度响应模方为式( 2-2) : (2-2) 式中 c w 为有效通带截止频率,表示与通带波纹有关的参量, 值越大通 带不动愈大, N表示滤波器阶数。 )(xCN是 N阶切比雪夫多项式,定义为(2-3)所示: (2-3) (三)切比雪夫低通滤波器的频域特性如式(2-4):、 间振荡和在时,)( 2 2 1 1 101jwww c 增大,下降加速)单调下降(时,)(Njwww c 2 2 控制了通带衰减)( 2 2 1 1 3 c jw j 1 c N=2 N=3 N=7 )/(1 1

12、 )j ( 22 2 cN C H 1)(harccoscosh 1 )arccos(cos )( xxN xxN xCN ss lg20A 6 2 2 2 1 1 0 10 jN jN 为偶数时, 为奇数时, (2-4) (四) 切比雪夫低通滤波器的设计步骤 )(arccos )110 1 (arccos 3 1102 1 1 .0 1.0 p s A A pcCp w w h h N N www s p 由通带、阻带指标确定)( 由通带衰减确定)( 确定由通带截频)( (2-5) 2双线性变换法 (一)基本思想 将非带限的模拟滤波器映射为最高频率的带限模拟滤波器, 其方法为式( 2-6 7

13、) : zss (2-6) TT ww Tw wT T w /,/, ) 2 arctan( 2 (2-7) (二) s 域到 z 域的映射关系如式( 2-810) : 2 tan 2 T w(2-8) 式中, T仍是采样间隔。 (2-9) ) 2 cos( ) 2 sin( 2 j) 2 tan( 2 jj TT 2 j 2 j 2 j 2 j ee ee2 T j j e1 e12 T 7 (2-10) (三)稳定性分析 令jws,则有(2-11) 域单位圆外域右半平面映射到 域单位圆上域虚轴映射到 域单位圆内域左半平面映射到 zSz zSz zSz 1,0)3( 1,0)2( 1,0)

14、1( 因此,稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。 (四)双线性变换法优缺点 双线性变换法与脉冲响应不变法相比,其主要的优点是避免了频率响应的混叠 现象。缺点是幅度响应不是常数时会产生幅度失真。 (五) 和的关系如下图 2-3 所示: 2 tan 2 T w 图 2-3 模拟频率 与数字频率 的对应关系 由图 2-3 看出,在零频率附近, 模拟角频率 与数字频率 之间的变换关系 接近于线性关系;但当 进一步增加时, 增长得越来越慢,最后当时, s T s T z z z T s T j ezjws 2 2 1 12 e1 e12 j 1 1 j j 、 1 1 1 12

15、 z z T s sT sT z 2 2 22 22 )/2( )/2( T T z /T /T 8 终止在折叠频率 =处,因而双线性变换就不会出现由于高频部分超过折叠 频率而混淆到低频部分去的现象,从而消除了频率混叠现象。 但是双线性变换的这个特点是靠频率的严重非线性关系而得到的,如式 (2-10)及图 2-3 所示。由于这种频率之间的非线性变换关系,就产生了新的问 题。首先,一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤 波器,不再保持原有的线性相位了;其次,这种非线性关系要求模拟滤波器的幅 频响应必须是分段常数型的, 即某一频率段的幅频响应近似等于某一常数(这正 是一般典型

16、的低通、高通、带通、带阻型滤波器的响应特性),不然变换所产生 的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变,如图2-4 所 示。 o o o )j ( a H )(e j H o o o )(earg j H )j (argaH 图 2-4 双线性变换法幅度和相位特性的非线性映射 对于分段常数的滤波器, 双线性变换后, 仍得到幅频特性为分段常数的滤波 器,但是各个分段边缘的临界频率点产生了畸变,这种频率的畸变, 可以通过频 率的预畸来加以校正。 也就是将临界模拟频率事先加以畸变,然后经变换后正好 映射到所需要的数字频率上。 9 三 IIR 数字滤波器设计过程 1. 设计步骤 根据以

17、上 IIR 数字滤波器设计方法,下面运用双线性变换法基于MATLAB 设计一个 IIR 带通滤波器。 (1) 确定性能指标 在设计带通滤波器之前, 首先根据工程实际的需要确定滤波器的技术指 标: 通带截止频率 fp1=1200, fp2=3000; 阻带截止频率 fs1=1000; fs2=3200; ; 阻带最小衰减 As=100dB 和通带最大衰减 Ap=1dB; (2) 把频率转化为数字角频率 wp1=2*pi*fp1*T; wp2=2*pi*fp2*T; ws1=2*pi*fs1*T; ws2=2*pi*fs2*T; (3) 频率预畸变 用=2/T*tan(w/2)对带通数字滤波器H(

18、z) 的数字边界频率预畸变, 得 到带通模拟滤波器H(s) 的边界频率主要是通带截止频率Wp1,Wp2; 阻带截止 频率 Ws1 ,Ws2的转换。抽样频率fs=10KHz。 通带截止频率 Wp1=(2/T)*tan(wp1/2); Wp2=(2/T)*tan(wp2/2); 阻带截止频率 Ws1=(2/T)*tan(ws1/2); Ws2=(2/T)*tan(ws2/2) (4) 模拟带通性能指标转换成模拟低通性能指标 WP=Wp1,Wp2;W0=sqrt(Wp1*Wp2); B=Wp2-Wp1; WS=Ws1,Ws2; (5) 模拟低通滤波器的构造 借助切比雪夫(Chebyshev) 滤波器

19、得到模拟低通滤波器的传输函数 Ha(s) 。 N,Wc=cheb1ord(WP,WS,Ap,As,s); %求阶数和边缘频率 z0,p0,k0=cheb1ap(N,Ap); %求极点,零点和增益 num=k0*real(poly(z0); %模拟低通滤波器系统函数的分子多项式 den=real(poly(p0); (6) 模拟低通滤波器转换成模拟带通滤波器 调用 lp2bp 函数将模拟低通滤波器转化为模拟带通滤波器。 numt,dent=lp2bp(num,den,W0,B); % 模拟带通滤波器系统函数的分子多项式和分母多项式 10 (7) 模拟带通滤波器转换成数字带通滤波器 利用双线性变换

20、法将模拟带通滤波器Ha(s) 转换成数字带通滤波器 H(z) 。 numd,dend=bilinear(numt,dent,fs); 2. 音频信号部分程序 (1)写入声音信号 y,fs,nbits=wavread (1.wav); sound(y,fs,nbits); %回放语音信号 N= length (y) ; %求出语音信号长度 Y=fft(y,N); %傅里叶变换 (2)产生噪声并加到声音中 noise=0.01*randn(N,1); %随机函数产生噪声 S=y+noise; %语音信号加入噪声 sound(S); Si=fft(S); %滤波前傅里叶变换 y1=filter(nu

21、md,dend,S); sound(y1); y2=fft(y1); 3. 程序流程图 首先确定性能指标, 把频率转化为数字角频率, 进而在进行频率预畸变, 用 =2/T*tan(w/2)对带通数字滤波器H(z) 的数字边界频率预畸变 , 得到带通模拟 滤波器 H(s) 的边界频率主要是通带截止频率Wp1,Wp2; 阻带截止频率Ws1 ,Ws2 的转换。抽样频率fs=10KHz。 上述准备工作做好之后,就先把模拟带通性能指标转换成模拟低通性能指 标,然后设计模拟低通滤波器,借助切比雪夫(Chebyshev) 滤波器得到模拟低通 滤波器的传输函数Ha(s) 。然后调用 lp2bp 函数将模拟低通

22、滤波器转化为模拟带 通滤波器。 最后利用双线性变换法将模拟带通滤波器Ha(s) 转换成数字带通滤波 器 H(z) 。 ;流程图如下图3-1 所示: 11 图 3-1 IIR带通滤波器去噪流程 . 仿真结果 源程序设计了模拟低通滤波器、模拟带通滤波器与数字带通滤波器,对数 字带通滤波器的性能仿真如图3-2: 00.10.20.30.40.50.60.70.80.91 -100 -50 0 w(rad) |H ( z ) | 幅频特性曲线(db) 00.10.20.30.40.50.60.70.80.91 -5 0 5 w(rad) H ( z ) 相频特性曲线 图 3-2 滤波器的辐频相频特性曲

23、线 开始 在 Window下录制语音将语音格式改为wav 对语音信号进行频谱分析,画出时域和频域波形图 加入噪声 画出其频率响应 用 IIR 滤波器对语音信号进行滤波 画出语音信号滤波前后的波形并进行比较 结束 用切比雪夫设计IIR 带通滤波器 12 为了实现 IIR 数字带通滤波器的应用,程序中加入了有噪声的音频信号, 通过对其滤波处理, 来显示数字带通滤波器的功能,下面显示未加入噪声如图 3-3 所示频谱及加入噪声之后和加噪后经滤除噪声的频谱和波形如图3-4 所 示: 0123456 x 10 4 -1 -0.5 0 0.5 1 原 始 信 号 波 形 0123456 x 10 4 0 5

24、00 1000 1500 2000 原 始 信 号 频 谱 图 3-3 语音信号的波形及频谱 0246 x 10 4 0 500 1000 1500 2000 滤波前信号频谱 0246 x 10 4 0 5 10 滤波后信号频谱 0246 x 10 4 -1 -0.5 0 0.5 1 滤波前信号波形 0246 x 10 4 -0.01 -0.005 0 0.005 0.01 滤波后信号波形 图 3-4 滤波前后的频谱与噪声 13 总结 在现代通信系统中, 由于信号中经常混有各种噪声和干扰,所以很多信号分 析都是基于滤波器而进行的, 而数字滤波器一定运算关系改变输入输出信号所含 频率成分的相对比

25、例或者滤出某些频率成分的器件,具有处理精度高、 稳定、灵 活、 不从在阻抗匹配问题等优点, 可以实现模拟滤波器无法实现的特殊滤波功能。 本次试验的主要利用MATLAB 软件实现基于语音信号去噪处理的IIR 滤波器的设 计。此次过程中主要有: 一、根据给定的指标设计得到IIR 数字带通滤波器; 二、 采集语音信号并对其进行时域、频域分析; 三、给采集好的信号叠加单频正弦噪 声;四、将叠加噪声之后的语音信号通过IIR 带通滤波器得到去噪语音信号,并 将去噪前后的语音信号进行时域、频域的分析对比,故而可以的到实验结果。 最后,虽然基于语音信号去噪处理的IIR 滤波器的课程设计得以完成, 但是也存 在

26、一些问题, 如选用不同的语音信号时, 经滤波后的不到与原始信号类似的信号 频谱,有事甚至什么语音信号都没有。因此,在以后的学习中要多加注意,并尽 可能找到解决的方法。 14 致谢 通过基于语音信号去噪处理的IIR 滤波器的设计我学习掌握了好多知识, 首 先是对 MATLAB 有了一个全新的认识, 其次是对 MATLAB 的更多操作和命令的使用 有了更高的掌握,最重要的是对MATLAB 处理数字信号处理的相关能力有了更高 的飞跃,就对 MATLAB 中涉及到数字信号处理的相关命令来说,通过这次课程设 计的亲身操作和实践,学习掌握了较多的原来不知道的或不熟悉的命令。 在此次基于语音信号去噪处理的I

27、IR 带通滤波器设计仿真过程中, 滤波器的 设计对于这次课设的我来说完全是实践中的新知识,虽然网上有很多关于IIR 带通滤波器的完整设计, 能有借鉴的模板和方法固然很令人放松,但自己觉得这 样做有点浪费青春, 很没意思。 因为我借鉴过来后, 自己对基于语音信号去噪处 理的 IIR 带通滤波器设计既没有记忆, 更谈不上理解。所以一个人暗暗下定决心, 要把这块硬骨头啃下来。 因此图书馆便是不可少去的地方,翻阅多本参考书, 进 行多重比较, 才在心中有了设计的基本概念。在多次爬山涉水中, 才明白怎样进 行完整的设计,怎样用MATLAB 实现相关的仿真。 最后, 非常感谢陈海燕老师和同学在这次基于语音

28、信号去噪处理的IIR 滤波 器课程设计中对我的帮助, 因为你们的不懈努力和耐心教导,才使在这次课程设 计中收获颇丰。 15 参考文献 1 陈后金 . 数字信号处理 M. 北京:高等教育出版社 ,2005.1 2 程佩青 . 数字信号处理 M. 北京:清华大学出版社 ,2007.2 3 从玉良 . 数字信号处理原理及其MATLAB 实现 M. 北京:电子工业出版 社.2009.7 4 胡广书 . 数字信号处理理论、 算法与实现 M. 北京:清华大学出版社 .2003.8 5 丁玉美,高西全 . 数字信号处理 M. 西安电子科技大学出版社.2001 6 刘泉. 数字信号处理原理与实现 M. 电子工

29、业出版社(第2 版).2009 7 陈怀琛 . 数字信号处理教程 MATLAB 释义与实现 M. 电子工业出版社 .2004 8 张圣勤 . MATLAB7.0实用教程 M. 北京: 机械工业出版社 .2006. 16 附录: % 技术指标 fp1=1200; fp2=3000; fs1=1000; fs2=3200; Ap=1; As=100; fs=10000; T=1/fs; wp1=2*pi*fp1*T; wp2=2*pi*fp2*T; ws1=2*pi*fs1*T; ws2=2*pi*fs2*T; % 带通到低通的频率转换 Wp1=(2/T)*tan(wp1/2); Wp2=(2/T

30、)*tan(wp2/2); WP=Wp1,Wp2; %模拟滤波器的通带截止频率 Ws1=(2/T)*tan(ws1/2); Ws2=(2/T)*tan(ws2/2); WS=Ws1,Ws2; %模拟滤波器的阻带截止频率 B=Wp2-Wp1; %带通滤波器的通带宽度 W0=sqrt(Wp1*Wp2); %带通滤波器的中心频率 % 切比雪夫模拟低通原型滤波器的设计 N,Wc=cheb1ord(WP,WS,Ap,As,s); %求阶数和边缘频率 z0,p0,k0=cheb1ap(N,Ap); %求极点,零点和增益 num=k0*real(poly(z0); %模拟低通滤波器系统函数的 分子多项式 d

31、en=real(poly(p0); %模拟低通滤波器系统函数的分母 多项式 % 模拟低通原型滤波器转换为模拟带通滤波器 numt,dent=lp2bp(num,den,W0,B); %模拟带通滤波器系统函数的 分子多项式和分母多项式 numd,dend=bilinear(numt,dent,fs); %双线性变换 , 有模拟滤波器 转为数字滤波器 h,w=freqz(numd,dend); %数字带通滤波器的幅频特性 dbH=20*log10(abs(h)+eps)/max(abs(h); 17 % 以分贝表示带通滤波器的幅频特性 figure(1); subplot(2,1,1);plot(

32、w/pi,db(h); grid on; xlabel(w(rad);ylabel(|H(z)|); axis(0,1,-100,25); title(幅频特性曲线 (db); subplot(2,1,2);plot(w/pi,angle(h);grid on; xlabel(w(rad); ylabel(H(z); axis(0,1,-5,5); title(相频特性曲线 ); y,fs,nbits=wavread (公鸡打鸣 .wav); sound(y,fs,nbits); %回放语音信号 N= length (y) ; %求出语音信号长度 Y=fft(y,N); %傅里叶变换 figu

33、re(2); subplot(2,1,1);plot(y);title(原始信号波形 ); subplot(2,1,2);plot(abs(Y); title(原始信号频谱 ); t=0:1/fs :1; noise=0.5*sin(2000*pi*t); %随机函数产生噪声 noise=noise(:,1); S=y+noise; %语音信号加入噪声 sound(S); Si=fft(S); %滤波前傅里叶变换 y1=filter(numd,dend,S); sound(y1); y2=fft(y1); %滤波后傅里叶变换 figure(3); subplot(2,2,1); plot(abs(Si),g); title(滤波前信号频谱 ); subplot(2,2,2); 18 plot(abs(y2),r); title(滤波后信号频谱 ); subplot(2,2,3); plot(S); title(滤波前信号波形 ); subplot(2,2,4); plot(y1); title(滤波后信号波形 );

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

当前位置:首页 > 其他


经营许可证编号:宁ICP备18001539号-1