数字信号处理-多采样率数字信号处理在数字语音系统中的应用

1、目录一、课程设计的性质与目的1二、课程设计题目11. 设计目的12. 设计要求13. 设计步骤2三、课程设计要求2四、设计进度安排2五、设计原理31. 巴特沃斯滤波器32. 采样定理43. 椭圆滤波器44. 抽取与内插5六、实验步骤及效果图61. 信源的时域和频域分析62. 对原始信号进行滤波73. 对滤波后的信号采样84. 椭圆滤波器滤波95. 对语音信号进行抽取和内插处理106. 语音信号的恢复117. 信号波形图对比12七、问题及解决办法141. 信源的时域和频域分析 原代码142. 对原始信号进行滤波 原代码153. 椭圆滤波器滤波 原代码16八、心得体会174一、课程设计的性质与目的

2、数字信号处理课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。通过该课程的课程设计实践，使学生对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解；巩固和运用在数字信号处理课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。2、 课程设计题目多采样率数字信号处

3、理在数字语音系统中的应用1、 设计目的学习多采样率数字信号处理原理，采用整数因子抽取与整数因子内插来解决数字语音系统中的信号采样过程中存在的问题，并用MATLAB编程实现，加深对多采样率数字信号处理的理解。2、 设计要求编程实现数字语音系统中改进的发送端A/D转换方案；编程实现数字语音系统中改进的接收端D/A转换器方案；给出各个设计环节信号的时域图及频域图。前置取样滤波A/D变换器数字信号处理器D/A变换器模拟滤波器模拟xa(t)模拟ya(t)采样量化转换为模拟电平DSP3、设计步骤(1) 设计模拟预滤波器，其过渡带为4f12 kHz。(2) 以16 kHz的采样率进行语音采样。(3) 设计数

4、字滤波器滤掉采样后频谱在412 kHz的频带中发生的混叠。(4) 按抽取因子D=2进行抽取，降低采样率，使得数据量不增加。(5) 按内插因子I=2进行内插，将采样率提高2倍。(6) 设计模拟恢复低通滤波器恢复的语音信号。3、 课程设计要求1、在一周内学生须上机16小时以上，程序调试完后，须由指导老师在机器上检查运行结果，经教师认可后的源程序可通过打印机输出，并请教师在程序清单上签字。2、课程设计报告内容和格式：设计题目，设计的详细步骤，设计过程中的结果、图形等，设计总结。3、每组每人必须独立完成，成绩的考核按设计结果、答辩成绩及课程设计报告来综合评定。成绩分为优、良、中、及格、不及格五级分评定

5、。4、指导教师：彭祯 曹燚 罗雅博。4、 设计进度安排通信工程1301/02:周周一上午，E-,布置课程设计任务，学生上机，按任务要求进行课程设计；分组选题；周周二上午，E-,学生上机，按任务要求进行课程设计；分组选题；周周三上午，E-,学生上机，按任务要求进行课程设计；分组选题；周周四上午，E-,学生任务完成，答辩并提交课程设计报告。5、 设计原理1、巴特沃斯滤波器设计巴特沃斯低通滤波器的方法：（1） 计算归一化频率；（2） 根据设计要求按照：计算巴特沃斯滤波器的参数C和阶次N，当ap=3db时，C=1；（3） 利用N查表获得归一化巴特沃斯低通原型滤波器的系统函数H（p）；（4） 令H（p）

6、中的得到截止频率为wp的巴特沃斯低通滤波器。设计巴特沃斯低通滤波器的指标：（1） wp-通带边界频率的归一化值； （2） ws-阻带边界频率的归一化值；（3） Rp-通带最大衰减；（4） Rs-阻带最小衰减。用MATLAB函数设计巴特沃斯滤波器：（1） Z,P,K=buttap(N);计算N阶巴特沃斯归一化模拟低通原型滤波器系统函数的零、极点增益因子；（2） N,wc=buttord(wp,ws,Rp,As)；（3） N,wc=buttord(wp,ws,Rp,As)；（4） B,A=butter(N,wc,ftype)；（5） B,A=butter(N,wc,ftype,s)。2、 采样定理

7、模拟信号经过A/D变换转换为数字信号的过程称之为采样，信号采样后其频谱产生了周期延拓，每隔一个采样频率ws，重复出现一次。为了保证采样后的信号能真实地保留原始模拟信号的信息，采样信号的频率必须至少为原信号中成分的2倍，这是采样的基本法则，成为采样定理。需要注意的是，在对信号进行采样时，满足采样定理只能保证不发生频率混叠，对信号的频谱做逆傅立叶变换时，可以完全变换为时域采样信号，而不能保证此时的采样信号能 真实反映原信号。一般频谱函数是复函数，相加应该是复数相加，只有当信号最高频率不超过Fs/2时，才不会产生频率混叠现象。3、 椭圆滤波器椭圆滤波器在通带和阻带内都具有等波纹幅频响应特性，是可以获

8、得对理想滤波器幅频响应的最好逼近，是一种性价比最高的滤波器。用MATLAB函数设计椭圆滤波器：（1） z,p,k=ellipap(N,Rp,As)；（2） N,wpo=ellipord(wp,ws,Rp,As)；（3） N,wpo=ellipord(wp,ws,Rp,As，s)；（4） B,A=ellip(N,Rp,wpo,ftype)；（5） B,A=ellip(N,Rp,wpo,ftype,s)；4、抽取与内插6、 实验步骤及效果图1、信源的时域和频域分析x,fs,bits=wavread('C:UserswlzDesktop2.wav'); %调用音频sound(x,fs

9、,bits); %声音x=x(:,1); %单声道Fs=fs %频率N=length(x); %音乐的长度t=0:N-1/fs;fk=0:fs/N:fs*(N-1)/N; y1=fft(x); %fftfigure(1)subplot(2,1,1)plot(t,x);title('原始语音信号时域波形');subplot(2,1,2)plot(fk,abs(y1);title('原始语音信号频谱')用wavread函数调用32khz的音频文件，用sound函数播放音乐，x=x(:,1)将其调成单声道音频，将原始信号的波形和频谱用matlab显示。2、对原始信号进

10、行滤波wp=2*4000/fs; %书P161页，巴特沃斯滤波器 过渡带为4KHz-12KHzws=2*12000/fs;Rp=1;Rs=40;N,wn=buttord(wp,ws,Rp,Rs); %wp-通带边界频率的归一化值 ws-阻带边界频率的归一化值 Rp-通带最大衰减 Rs-阻带最小衰减B,A=butter(N,wn);K=0:511;fk1=0:fs/512:fs*(511/512); %均匀取512个点Hk=freqz(B,A,512);%计算频率响应figure(2);subplot(111)plot(fk1',20*log10(abs(Hk);grid onxlabe

11、l('频率/Hz');ylabel('频率响应幅度');title('Butterworth')y=filter(B,A,x);Yk=fft(y);figure(3)subplot(2,1,1)plot(t,y);title('滤波后信号时域波形');subplot(2,1,2)plot(fk,abs(Yk);title('滤波后语音信号频谱');利用巴特沃斯滤波器设计一个过渡带为4KHz-12KHz的滤波器，wp=2*4000/fs;ws=2*12000/fs;Rp=1;Rs=40;滤波后的信号波形和频谱如图所示

12、。3、对滤波后的信号采样fs=16000; %以16kHz进行采样N=length(x);n=0:length(x)-1;T=1/Fs;Tp=(N-1)*T;ta=0:(Fs/fs)*T:Tp;xat=x(1:Fs/fs:N); Na=length(xat);y3=fft(xat);fka=0:fs/Na:fs*(Na-1)/Na;sound(xat,fs,bits)figure(4)subplot(2,1,1)plot(ta,xat)grid on xlabel('t/s');ylabel('y2(t)');title('采样后语音信号的波形'

13、)grid on;subplot(2,1,2)plot(fka,abs(y3)/max(abs(y3)title('采样后语音信号的频率')xlabel('f/Hz');ylabel('幅度');以16KHz对滤波后的信号进行采样，采样结果如图所示。4、椭圆滤波器滤波N,wn=ellipord(wp,ws,Rp,Rs);%计数字滤波器滤掉采样后频谱在412 kHz的频带中发生的混叠。B,A=ellip(N,Rp,Rs,wn); %椭圆滤波器h,w=freqz(B,A);figure(4)subplot(1,1,1)grid on;plot(w*f

14、s/(2*pi),20*log10(abs(h)/max(abs(h); title('数字低通滤波器的幅度响应');xlabel('频率/Hz');ylabel('频率响应幅度');title('ellipor')xx=filter(B,A,xat); %滤波sound(xx,fs,bits);figure(5)subplot(2,1,1)stem(ta,xx,'.');axis(1.11,1.112,-1,1)title('椭圆滤波后的信号时域波形')yy1=fft(xx);subplot(2,

15、1,2)plot(abs(yy1);title('椭圆滤波后的信号频谱')用椭圆滤波器使采样后频谱在412 kHz的频带中发生的混叠。5、对语音信号进行抽取和内插处理y4=decimate(xx,2); %抽取，对时间序列进行整数倍采样处理，使得长度降低，采样率减低为原来的1/2y5=interp(y4,2); %内插，对时间序列进行整数倍插值，使得时间序列曲线更光滑，插值后使得序列变为2倍y41=fft(y4);y51=fft(y5);sound(y4,fs,bits);figure(6)subplot(2,2,1);dt=0:4*T:Tpstem(dt,y4,'.&

16、#39;)axis(1.11,1.112,-1,1)grid on xlabel('2倍抽取后的信号');subplot(2,2,2);plot(abs(y41)grid onxlabel('2倍抽取后的语音信号频谱');subplot(2,2,3);dt1=0:2*T:Tpstem(dt1,y5,'.')axis(1.11,1.112,-1,1)grid onxlabel('2倍插值后的信号');subplot(2,2,4);plot(abs(y51)grid onxlabel('2倍插值后语音信号频谱');(7

17、) 两倍抽取和两倍插值后的信号波形及频谱。过程中使用decimate函数对信号进行抽取，用interp函数对信号进行内插。按抽取因子D=2进行抽取，降低采样率，使得数据量不增加，按内插因子I=2进行内插，将采样率提高2倍。6、语音信号的恢复figure(7)subplot(2,2,1)plot(fk,x);axis(14800,15000,-1,1)title('语音信号原有32k采样波形');subplot(2,2,2);plot(abs(y1);title('语音信号原有32k采样波形频谱');n=1:length(xat);n2=1:2:length(xa

18、t);yi=interp1(n,xat,n2);sound(yi);subplot(2,2,3)plot(n2,yi);axis(14800,15000,-1,1)title('内差恢复波形');y2=fft(yi);subplot(2,2,4);plot(abs(y2);title('语音信号内差恢复波形频谱');使用interp1函数对信号进行恢复。7、 信号波形图对比原始语音信号巴特沃斯滤波器滤波后信号采样后语音信号由于使用的音频杂音较少，所以在滤波后的信号与原信号相比，波形变化不大，以16KHz频率采样使得在原有的32KHz的基础上采样点减少一半，采样后

19、的语音信号波形要比原信号波形稀疏。椭圆滤波后的信号 2倍抽取和2倍采样后的信号恢复波形在2倍抽取后，信号的采样点较椭圆滤波器滤波后减少，2倍内插后，采样点增多，内插恢复波形大致与原始信号相同，在播放中已基本和原信号相同。7、 问题及解决办法在实验中也并非是一帆风顺的，在过程中我遇到了很多问题，幸好有彭老师的耐心指导，才使得我将实验做出。1. 信源的时域和频域分析 原代码x,fs,bits=wavread('C:UserswlzDesktop2.wav'); %调用音频Fs=fsN=length(x);fk=0:fs/N:fs*(N-1)/N;y1=fft(x); %fftfig

20、ure(1)subplot(2,1,1)plot(x);title('原始语音信号时域波形');subplot(2,1,2)plot(abs(y1);title('原始语音信号频谱')在此代码中，虽然可以显示出信号的波形和频谱，但并未设置时间t,导致横坐标出现了问题，无法清晰分辨出之后各步骤执行完成后的结果是否正确，在彭老师的指导后，在程序中加入t=0:N-1/fs,plot(t,x)解决了该问题,同时用x=x(:,1)将信号变为单声道，使信号一目了然，没有重叠。2. 对原始信号进行滤波 原代码 wp=2*4000/fs; %设计模拟滤波器，过渡带为4<=

21、f<=12kHzws=2*12000/fs;Rp=2;Rs=30;N,wn=buttord(wp,ws,Rp,Rs,'s');B,A=butter(N,wn,'s');K=0:511;fk=0:fs/512:fs*(511/512);Hk=freqs(B,A,512);figure(2);subplot(111)plot(fk,20*log10(abs(Hk);grid onxlabel('频率/Hz');ylabel('频率响应幅度');title('Butterworth')y=filter(B,A,x)

22、;Yk=fft(y);figure(3)subplot(2,1,1)plot(t,y);title('滤波后信号时域波形');subplot(2,1,2)plot(fk,abs(Yk);title('滤波后语音信号频谱');在源代码中，图形并没有达到4KHz-12KHz过渡带预期的效果，在老师的指导下，我们将第五行改为N,wn=buttord(wp,ws,Rp,Rs)，将模拟滤波器改为数字滤波器，将Rp设为1，将Rs设为40，将Hk=freqs(B,A,512);改为Hk=freqz(B,A,512);最后得出的结果大致为正确结果。3.椭圆滤波器滤波 原代码N,

23、wn=ellipord(wp,ws,Rp,Rs); %椭圆滤波器,wp通带边界频率，ws阻带边界频率，Rp通带最大波纹，Rs阻带最小衰减B,A=ellip(N,Rp,Rs,wn);h,w=freqz(B,A);figure(4)subplot(1,1,1)grid on;plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h);title('数字低通滤波器的幅度响应');xlabel('频率/Hz');ylabel('频率响应幅度');title('ellipor')xx=filter(B,A,xat

24、);figure(5)subplot(2,1,1)plot(xx);title('滤波后的信号时域波形')yy1=fft(xx);subplot(2,1,2)plot(abs(yy1);title('滤波后的信号频谱')在椭圆滤波器的使用中，目的是滤掉采样后频谱在412 kHz的频带中发生的混叠，本段问题是显示的图形中横坐标取值范围太大，导致看不清图像，老师建议我们将横坐标范围缩小，以看清图形，我们在源代码的基础上，加入stem(ta,xx,'.')、axis(1.11,1.112,-1,1)，在之后的抽取、内插和恢复中，也使用该函数显示图形，使

25、图形一目了然，使横坐标保持一致。8、 心得体会 附录：x,fs,bits=wavread('C:UserswlzDesktop2.wav'); %调用音频sound(x,fs,bits); %声音x=x(:,1); %单声道Fs=fs %频率N=length(x); %音乐的长度t=0:N-1/fs;fk=0:fs/N:fs*(N-1)/N; y1=fft(x); %fftfigure(1)subplot(2,1,1)plot(t,x);title('原始语音信号时域波形');subplot(2,1,2)plot(fk,abs(y1);title('原始

26、语音信号频谱')wp=2*4000/fs; %书P161页，巴特沃斯滤波器 过渡带为4KHz-12KHzws=2*12000/fs;Rp=1;Rs=40;N,wn=buttord(wp,ws,Rp,Rs); %wp-通带边界频率的归一化值 ws-阻带边界频率的归一化值 Rp-通带最大衰减 Rs-阻带最小衰减B,A=butter(N,wn);K=0:511;fk1=0:fs/512:fs*(511/512); %均匀取512个点Hk=freqz(B,A,512);%计算频率响应figure(2);subplot(111)plot(fk1',20*log10(abs(Hk);gri

27、d onxlabel('频率/Hz');ylabel('频率响应幅度');title('Butterworth')y=filter(B,A,x);Yk=fft(y);figure(3)subplot(2,1,1)plot(t,y);title('滤波后信号时域波形');subplot(2,1,2)plot(fk,abs(Yk);title('滤波后语音信号频谱');fs=16000; %以16kHz进行采样N=length(x);n=0:length(x)-1;T=1/Fs;Tp=(N-1)*T;ta=0:(Fs/

28、fs)*T:Tp;xat=x(1:Fs/fs:N); Na=length(xat);y3=fft(xat);fka=0:fs/Na:fs*(Na-1)/Na;sound(xat,fs,bits)figure(4)subplot(2,1,1)plot(ta,xat)grid on xlabel('t/s');ylabel('y2(t)');title('采样后语音信号的波形')grid on;subplot(2,1,2)plot(fka,abs(y3)/max(abs(y3)title('采样后语音信号的频率')xlabel(

29、9;f/Hz');ylabel('幅度');N,wn=ellipord(wp,ws,Rp,Rs);%计数字滤波器滤掉采样后频谱在412 kHz的频带中发生的混叠。B,A=ellip(N,Rp,Rs,wn); %椭圆滤波器h,w=freqz(B,A);figure(4)subplot(1,1,1)grid on;plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h); title('数字低通滤波器的幅度响应');xlabel('频率/Hz');ylabel('频率响应幅度');title('ellipor')xx=filter(B,A,xat); %滤波sound(xx,fs,bits);figure(5)subplot(2,1,1)stem(ta,xx,'.');axi