数字滤波是数字信号处理的重要内容，数字滤波器可分为IIR和FIR两大类。对于IIR数字滤波器的设计，需要借助模拟原型滤波器，再将模拟滤波器转化为数字滤波器，文中采用的设计方法是脉冲响应不变法、双向性变换法和完全函数设计法；对于FIR数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。根据IIR滤波器和FIR滤波器的特点，在MATLAB坏境下分别用双线性变换法设计IIR和用窗函数设计FIR数字滤波器，并对采集的语音信号进行分析，最后给出了IIR和FIR对语音滤波的效果。
1.2 课题要求
1. 掌握数字信号处理的基本概念，基本理论和基本方法。 2. 熟悉离散信号和系统的时域特性。 3. 掌握序列快速傅里叶变换方法。
4. 学会MATLAB的使用，掌握MATLAB的程序设计方法。 5. 掌握利用MATLAB对语音信号进行频谱分析。 6. 掌握滤波器的网络结构。
第 1页 共39页
语音信号的数字滤波处理(九)
2 课程设计预习与原理
2.1 卷积运算的演示
2.1.1 线性卷积
序列x1(n)=[2 0 1 2 5 7 0 5 0 2 0 3]，序列x2(n)=[2 0 0 1 0 1 1 0]。动态演示两个序列进行线性卷积x1(n)﹡x2(n)的翻转、移位、乘积、求和的过程。其中翻转采用fliplr，程序如下：
n=-7:18; M=17;
yn=zeros(1,19); figure(1) stem(yn); xlabel('n') ylabel('y(n)')
xn1=[2 0 1 2 5 7 0 5 0 2 0 3];
xm1=[zeros(1,7) xn1 zeros(1,7)];%为26个数字的矩阵 figure(2) stem(n,xm1) xlabel('m') ylabel('x1(m)') xn2=[2 0 0 1 0 1 1 0];
xm2=[fliplr(xn2) zeros(1,18)]; %移位，补零为26个数字的矩阵 figure(3) stem(n,xm2) xlabel('m') ylabel('x2(N-m)') title('n=0')
yn(1)=sum((xm1.*xm2)');%对xm1与xm2进行对应原素乘方之后进行数组转置，求和;即为求卷积
figure(4)
stem(yn) xlabel('n') ylabel('y(n)') title('n=N') for N=1:17
xm3=[zeros(1,N) fliplr(xn2) zeros(1,M)]; figure(); stem(n,xm3)
第 1页 共39页
语音信号的数字滤波处理(九)
xlabel('m') ylabel('x2(N-m)') title('n=N') M=M-1;
yn(N+1)=sum((xm1.*xm3)'); figure() stem(yn) xlabel('n') ylabel('y(n)') title('n=N') end
xm3=[zeros(1,18) fliplr(xn2)] figure() stem(xm3) xlabel('m') ylabel('x2(N-m)') title('n=N');
yn(19)=sum((xm1.*xm3)'); figure() stem(yn) xlabel('n') ylabel('y(n)')
2.1.2循环卷积
序列x1(n)=[2 0 1 2 5 7 0 5 0 2 0 3]，序列x2(n)=[2 0 0 1 0 1 1 0]，N=12。动态演示两个序列进行圆周卷积x1(n)⊙x2 (n)的翻转、移位、乘积、求和的过程。程序如下：
n=0:11;
yn=zeros(1,12); figure(1) stem(yn) xlabel('n')
ylabel('y(n)') %图1,12个0
xn1=[2 0 1 2 5 7 0 5 0 2 0 3]; figure(2) stem(n,xn1) xlabel('m')
ylabel('x1(m)')
xn2=[2 0 0 1 0 1 1 0];
xm2=[xn2 zeros(1,length(xn1)-length(xn2))]; figure(3) stem(n,xm2) xlabel('m')
第 1页 共39页
语音信号的数字滤波处理(九)
ylabel('x2(m)') title('n=0');
yn(1)=sum((xn1.*xm2)'); figure(4) stem(yn) xlabel('n') ylabel('y(n)') title('n=N')
for N=1:11
xm1=[fliplr(xn1(1:N)) fliplr(xn1(N+1:12))]; figure() stem(n,xm1) xlabel('m') ylabel('x1(N-m)') title('n=N')
yn(N)=sum((xm1.*xm2)'); figure() stem(n,yn) xlabel('n') ylabel('y(n)') title('n=N') end
figure()
xm1=fliplr(xn1); stem(n,xm1) xlabel('m') ylabel('x1(N-m)') title('n=N')
yn(12)=sum((xm1.*xm2)'); figure() stem(n,yn) xlabel('n') ylabel('y(n)') title('n=N')
2.1.3声音文件线性卷积
序列x1(n)=[2 0 0 1 0 1 1 0],读取一段声音数据，当循环卷积长度大于或等于两序列长度之和时，循环卷积等于线性卷积。因为直接用FFT进行1024点卷积大于两序列长
第 1页 共39页
语音信号的数字滤波处理(九)
度，所以可用线性卷积替代圆周卷积，其程序如下：
n=0:686; M=2053;
yn=zeros(1,687);
[xn1,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); xn1=xn1(:); xn1=xn1';
xn1=xn1(10000:12047); %xn1=xn1';
xm1=[zeros(1,7) xn1 zeros(1,7)]; xn2=[2 0 0 1 0 1 1 0];
xm2=[fliplr(xn2) zeros(1,2054)]; yn(1)=sum((xm1.*xm2)'); for N=1:685
xm3=[zeros(1,N) fliplr(xn2) zeros(1,M)]; M=M-1;
yn(N+1)=sum((xm1.*xm3)') end
xm3=[zeros(1,2054) fliplr(xn2)]; yn(687)=sum((xm1.*xm3)'); figure(1) stem(yn) xlabel('m') ylabel('y(n)') title('n=N')
线性卷积结果如图所示
第 1页 共39页
语音信号的数字滤波处理(九)
图1
由循环卷积定理可知：对于时域序列循环卷积，可先进行FFT变换，然后频率相乘，最后对结果进行IFFT变换，即可得到时域循环卷积结果。其程序如下：
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); y=y(10000:12047); Y=fft(y,1024);
xn2=[1 2.43 6.17 12.93 22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.93 6.17 2.43 1]; X=fft(xn2,1024); X1=rot90(X,3); Z=(X1'.*Y)'; z=ifft(Z,1024)*8 figure(2) stem(z)
axis([0,700,-1,0.6]);
结果如图所示：
第 1页 共39页
语音信号的数字滤波处理(九)
图2
2.1.4采样定理的演示
序列长度为28，先计算x(n)的512点FFT，得到其频谱函数X512(k)，再对X512(k)隔点抽取，得到X32(k)和X16(k)。再分别对其进行IFFT变换。观察其时域图。其程序如下：
fs=27; n=3; t=0:1/fs:1;
xa=n*exp(-n*sqrt(2)*pi)*sin(sqrt(2)*n*t); subplot(3,2,1) stem(xa)
title('(a) 原28点时域信号xa') Xk=fft(xa,512); subplot(3,2,2) stem(Xk)
title('(b) 512点FFT[x(t)]') X32k=Xk(1:16:512); subplot(3,2,4) stem(X32k)
title('(d) 频率采样频率为32时的频谱图') x32a=ifft(X32k); subplot(3,2,3)
第 1页 共39页
语音信号的数字滤波处理(九)
stem(x32a);
title('(c) 32点IFFT[X32(k)]得到的x32a(t)') X16k=X32k(1:2:32); subplot(3,2,6) stem(X16k)
title('(f) 频率采样频率为16时的频谱图') x16a=ifft(X16k); subplot(3,2,5) stem(x16a);
title('(e) 16点IFFT[X16(k)]得到的x16a(t)')
结果如图所示：
图3
由图可知：当进行32点IFFT时，无时域混叠失真；当进行16点IFFT时，有时域混叠失真。
第 1页 共39页
语音信号的数字滤波处理(九)
2.2 课程设计原理
2.2.1 IIR设计原理 (1) 切比雪夫滤波器原理
切比雪夫滤波器的幅频特性具有等波纹特性。它有两种形式：振幅特性在通带内是等波纹的，在阻带内是单调下降的切比雪夫I型；振幅特性在阻带内是等波纹的，在通带内是单调下降的切比雪夫II型。
(2) 双线性变换法[4]工作原理
双线性变换中数字域频率和模拟频率之间的非线性关系了它的应用范围，只有当非线性失真是允许的或能被忽略时，才能采用双线性变换法，通常低通、高通、带通和带阻等滤波器等具有分段恒定的频率特性，可以采用预畸变的方法来补偿频率畸变，因此可以采用双线性变换设计方法。
(3) 脉冲响应不变法[4]工作原理
冲激响应不变法遵循的准则是使数字滤波器的单位取样响应与参照的模拟滤波器的脉冲响应的取样值完全一样，即h(n)=ha(nT),其中T为取样周期。实际是由模拟滤波器转换成为数字滤波器，就是要建立模拟系统函数Ha（S）与数字系统函数H（z）之间的关系。脉冲响应不变法是从S平面映射到z平面，这种映射不是简单的代数映射，而是S平面的每一条宽为2π/T的横带重复地映射到整个z平面。
2.2.2 FIR设计原理
由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。但是这些特性是以牺牲线性相位频率特性为代价的，即用Butterworth、切比雪夫和椭圆法设计的数字滤波器逼近理想的滤波器的幅度频率特性，得到的滤波器往往是非线性的。在许多电子系统中，对幅度频率特性和线性相位特性都有较高的要求，所以IIR滤波器在这些系统中往往难以胜任。
第 1页 共39页
语音信号的数字滤波处理(九)
３ 设计程序的调试和运行结果
3.1声音信号的提取与加噪试和运行结果
3.1.1提取声音信号
x1=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav',10); x1=x1';x1=x1(1,:);
x2=[1 2.43 6.17 12.93 22.17 32.25 40.88 45.87 45.87 40.88 32.25 22.17 12.93 6.17 2.43 1];
subplot(7,1,1); stem(x1); ylabel('x1(n)'); title('x1(n)'); subplot(7,1,2); stem(x2); ylabel('x2(n)'); title('x2(n)'); y=conv(x1,x2);
subplot(7,1,3); stem(y); ylabel('y');
title('x1(n)与x2(n)的卷积'); N1=length(x1); N2=length(x2); N=N1+N2-1;
X1=fft(x1,N);X2=fft(x2,N); subplot(7,1,4); stem(X1); ylabel('X1');
title('x1(n)的N点DFT'); subplot(7,1,5); stem(X2); ylabel('X2');
title('x2(n)的N点DFT'); Y1=X1.*X2;
第 1页 共39页
语音信号的数字滤波处理(九)
subplot(7,1,6); stem(Y1); ylabel('Y1');
title('X1与X2相乘的结果'); Y2=ifft(Y1); subplot(7,1,7); stem(Y2); ylabel('Y2');
title('Y1的IDFT结果');
运行结果如图：
图4 原始信号波形 ３.１.２ 对声音信号加噪 MATLAB程序如下：
[y,fs,bits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); %x：语音数据；fs：采样频率；bits：采样点数 sound(y,fs,bits); %话音回放 n=length (y); %求出语音信号的长度 Y=fft(y,n); %傅里叶变换 subplot(2,1,1); plot(y);
title('原始信号波形'); subplot(2,1,2); plot(abs(Y));
第 1页 共39页
语音信号的数字滤波处理(九)
title('原始信号频谱'); L=length(y);
noise=0.5*randn(L,2); y_z=y+noise;
sound(y_z,fs,bits);
m=length (y_z); %求出被污染语音信号的长度 Y_Z=fft(y_z,n); %傅里叶变换 subplot(3,1,1); plot(y_z);
title('被污染信号波形'); subplot(3,1,2); plot(abs(Y_Z));
title('被污染的语音信号频谱');
运行结果如图：
图5 被污染的信号
3.2 切比雪夫低通滤波器程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavresd(‘filename’)读取语音信号数据。其程序如下： [y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋
-你把我灌醉.wav')
第 1页 共39页
语音信号的数字滤波处理(九) 图6 原始信号 如图所示为加入噪声后信号的时域图和频谱图。其主要程序如下：x1=x+zs'
图7 加入噪声后信号
如图所示为低通滤波器的幅频图。设计指标：wpz=1000Hz; wsz=1200Hz; rp=1; rs=100;其设计程序如下：
第 1页 共39页
语音信号的数字滤波处理(九) 图8 切比雪夫低通滤波器
如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：
yd=filter(az,bz,x1);
图9 滤波后信号
第 1页 共39页
语音信号的数字滤波处理(九)
3.3切比雪夫高通滤波器程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(‘’)读取语音信号数据。其程序如下：
[x,fs,nbits]=wavread(‘C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav’)
图10 原语音信号的时域图和频谱图
如图所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：
x1=x+zs1;
图11 被污染语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九) 如图所示为高通滤波器的幅频图。设计指标：wpz=500Hz; wsz=700Hz; rp=1; rs=100;其设计程序如下： wpz=2*700/fs;wsz=2*500/fs;
wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=1;rs=100; [N,wc]=cheb2ord(wp,ws,rp,rs,'s'); [b,a]=cheby2(N,rs,wc,'high','s'); [bz,az]=bilinear(b,a,fs); Hk=freqz(bz,az,16384,fs);
图12 高通滤波器的幅频图
如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：
yd=filter(az,bz,x1);
图13 滤波后语音信号的时域图和频谱图
3.4 切比雪夫带通程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(‘filename’)读取语音信号数据。其程序如下：
第 1页 共39页
语音信号的数字滤波处理(九)
[x,Fs,bits]=wavread('C:\\Users\\acer\\Desktop\\ca9e36e190eb780c00b59c4459a6d5.wav');
图14 原语音信号的时域图和频谱图
如图所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：
x1=x+zs1;
图15 被污染语音信号的时域图和频谱图
如图所示为带通滤波器的幅频图。设计指标：fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB;其设计程序如下：
wpz=[2*1200/fs,2*3000/fs]; wsz=[2*1000/fs,2*3200/fs]; wp=2*fs*tan(wpz*pi/2); ws=2*fs*tan(wsz*pi/2); rp=1;rs=100;
第 1页 共39页
语音信号的数字滤波处理(九)
[N,wc]=cheb2ord(wp,ws,rp,rs,'s'); [b,a]=cheby2(N,rs,wc,'s'); [bz,az]=bilinear(b,a,fs); Hk=freqz(bz,az,16384,fs);
图16 带通滤波器的幅频图
如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：
yd=filter(az,bz,x1);
图17 滤波后语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九) ３.５ hamming低通程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(‘filename’)读取语音信号数据。其程序如下：
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav');
图18 原语音信号的时域图和频谱图
如图所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：z=x+y;
第 1页 共39页
语音信号的数字滤波处理(九)
图19 被污染语音信号的时域图和频谱图
如图所示为低通滤波器的幅频图和相频图。设计指标：wpz=1000Hz; wsz=1200Hz; rs=100;其设计程序如下：
wpz=2*1000*pi/fs;wsz=2*1200*pi/fs;rs=100; Bt=wsz-wpz;
N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2); wc=(wpz+wpz)/2/pi;
hn=fir1(N-1,wc,hamming(N)); figure(3)
freqz(hn,1,32678);
第 1页 共39页
语音信号的数字滤波处理(九) 图20 低通滤波器的幅频图和相频图 如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下： hn=filter(hn,1,z);
图21 滤波后语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九)
３.６ hamming高通程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(‘filename’)读取语音信号数据。其程序如下：
[[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav');
图22 原语音信号的时域图和频谱图
如图所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：
z=x+y;
图23 被污染语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九) 如图所示为高通滤波器的幅频图和相频图。设计指标：wpzl=2800Hz; wszl=3200Hz; rs=100;其设计程序如下： wpz=2*3000*pi/fs;wsz=2*2800*pi/fs; rs=150; Bt=wpz-wsz;
N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2); wc=(wpz+wpz)/2/pi;
hn=fir1(N-1,wc,'high',hamming(N));
图24 高通滤波器的幅频图和相频图
如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：
hn=filter(Hn,1,z);
第 1页 共39页
语音信号的数字滤波处理(九) 图25 滤波后语音信号的时域图和频谱图 ３.７ hamming带通程序的调试和运行结果
如图所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(‘filename’)读取语音信号数据。其程序如下：
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav');
图26 原语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九) 如图所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下： t=0:0.000045:0.90855;
x1=0.1*sin(2*pi*10*t)+0.2*sin(2*pi*20*t)+0.07*sin(2*pi*6000*t); x=[rot90(x1,2),rot90(x1,2)]; x=x(1:32768); z=x+y;
图27 被污染语音信号的时域图和频谱图
如图所示为带通滤波器的幅频图和相频图。设计指标：wpz=4500Hz; wsz=4700Hz; wpzl=700Hz; wszl=500Hz; rp=20; rs=50;其设计程序如下：
wpzl=1200*pi/fs;wpzu=3000*pi/fs;wszl=1000*pi/fs;wszu=3200*pi/fs;rs=100; Bt=wpzl-wszl; N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2);
wc=[(wpzl+wszl)/2/pi,(wpzu+wszu)/2/pi]; Hn=fir1(N-1,wc,hamming(N));
第 1页 共39页
语音信号的数字滤波处理(九) 图28 带通滤波器的幅频图和相频图 如图所示为滤波后语音信号的时域图和频谱图。其滤波程序如下： hn=filter(hn,1,z);
图29 滤波后语音信号的时域图和频谱图
第 1页 共39页
语音信号的数字滤波处理(九)
４ 总结及心得体会
近一周的课程设计终于接近尾声，在这一周里，我感觉到了自己很多的不足，同时也有着诸多的辛酸与收获。
首先，这次的课程设计与以往的课程设计有着很大的不同之处，在时间上相对于以往两周的课设，这次的时间仅仅为一个星期，或者说4天，这个对我们来说时间上还是比较紧的，需要合理的安排分配时间才能够及时的完成；在分组上，此次的设计采用一人一组的形式，和以往一组2~3人的形式相比，这次的实习就缺乏了组员之间的分工与讨论，对个人解决问题的能力还是要求相对比较高的。总的来说，这次的实习对我来说是一项挑战吧，需要我迎难而上认真完成。
在课设的过程中让我感受很深的还是自己能力的不足，让自己觉得挺失望的。这次的实习主要是以数字信号处理的知识为理论基础，以MATLAB编程为工具完成对信号的处理。可是自我反省，觉得首先在数字信号的理论基础这一点上，那些理论知识基本上是学过之后便忘记了，不会付诸于应用，然后再课设的过程中需要自己去花费大量的时间学习，所以就会想要啥自己的理论知识足够扎实就好了。其次对于MATLAB的应用，觉得以前对于MATLAB的学习还是不够深入，基本上是依样画葫芦的设计一些较为基础简单的程序，不能灵活的运用于自我编程。可是现在的课设和以往的平时的学习就不一样了，需要自己灵活的运用所学的知识进行编程，不得不说让我觉得很烧脑，有时候一个程序调试了大半天都没有得到想要的结果，那种感觉确实很不好。
虽说课设烧脑那是肯定的，但是收获也不少。我觉得课设的目的之一就是以课程设计来实践完善平时理论知识的不足，同时锻炼动手操作能力，可以说课程设计和理论学习是一个相辅相成的过程。我也确实对MATLAB这样一个强大的工具有了更深刻的理解，比如其对声音信号的提取，加入噪声的过程以及滤除噪声的过程。若联系实际生活我想这与我们生活中的声音、图像等方面信息的处理也是有着很大的关联的，不得不说其实在大学里若认真去学，还是可以学到很多精彩的东西。我很清楚的知道自己是本科毕业以后将要就要面临找工作的人，通过这样一次课程设计来补充自己理论知识的不足，提高自己的动手操作能力，为自己以后的就业夯实基础，我觉得这是一个很有效的途径。虽然自己现在的能力与以后就业所需要的能力相比还存在着较大的差距，但是我
第 1页 共39页
语音信号的数字滤波处理(九)
也不可以因此而自暴自弃，不能让自己的能力不足成为我停滞不前的借口。学习是一个终身的的过程，虽然这次的课设已经结束，我想其在我以后学习生活道路上的影响必定是深远的。
第 1页 共39页
语音信号的数字滤波处理(九)
参考文献
[1] 刘卫国．MATLAB程序设计应用[M]．北京：高等教育出版社，2008． [2] 赵红怡．数字信号处理及MATLAB实现[M]．北京：化学工业出版社，2002． [3] 高西全，丁玉美．数字信号处理（第三版）[M]．陕西：西安电子科技大学出版社，2010．
[4] 何强，何英．MATLAB扩展编程[M]．北京：清华大学出版社，2002．
第 1页 共39页
语音信号的数字滤波处理(九)
附录A 切比雪夫滤波器完整程序
I切比雪夫低通滤波器
[x,Fs,bits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); x=x(:,1); figure(1); subplot(2,1,1); plot(x);
sound(x,Fs,bits);
title('语音信号时域波形图'); y=fft(x,3260); f=(Fs/1630)*[1:1630]; subplot(2,1,2);
plot(f(1:1630),abs(y(1:1630))); title('语音信号频谱图'); t=0:length(x)-1;
zs=0.05*cos(2*pi*10000*t/22050); zs0=0.05*cos(2*pi*10000*t/22050000); figure(2); subplot(2,1,1) plot(zs0)
title('噪声信号波形'); zs1=fft(zs,1200);
sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)
plot(f(1:600),abs(zs1(1:600))); title('噪声信号频谱'); %加噪声之后的语音 x1=x+zs';
%sound(x1,Fs,bits); 回放加入噪声后的语音 y1=fft(x1,1200); figure(3);
subplot(3,1,1);plot(x1); title('加入噪声后的信号波形'); subplot(3,1,2);
plot(f(1:600),abs(y1(1:600))); title('加入噪声后的信号频谱') %低通滤波器
wpz=2*4500/fs;wsz=2*4700/fs;
wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=2;rs=35; [N,wc]=cheb2ord(wp,ws,rp,rs,'s');
第 1页 共39页
语音信号的数字滤波处理(九)
[b,a]=cheby2(N,rs,wc,'s'); [bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(4);
plot(w,abs(h));
title('切比雪夫低通滤波器'); xlabel('频率（HZ）'); ylabel('耗损（dB）'); grid on;
yd=filter(az,bz,x1); figure(5);
subplot(2,1,1); plot(yd);
title('滤波后的时域波形图'); ydd=fft(yd,800); subplot(2,1,2);
plot(f(1:600),abs(ydd(1:600))); title('滤波后的频域波形图'); sound(yd,Fs,bits);
II切比雪夫高通滤波器
[x,Fs,bits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); x=x(:,1); figure(1); subplot(2,1,1); plot(x);
sound(x,Fs,bits);
title('语音信号时域波形图'); y=fft(x,3260); f=(Fs/1630)*[1:1630]; subplot(2,1,2);
plot(f(1:1630),abs(y(1:1630))); title('语音信号频谱图'); %被污染的语音 L=length(x); zs0=0.5*randn(L,2); figure(2); subplot(2,1,1) plot(zs0)
title('噪声信号波形'); zs1=fft(zs,1200);
第 1页 共39页
语音信号的数字滤波处理(九)
sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)
plot(f(1:600),abs(zs1(1:600))); title('噪声信号频谱'); %加噪声之后的语音 x1=x+zs';
%sound(x1,Fs,bits); %回放加入噪声后的语音 y1=fft(x1,1200); figure(3);
subplot(3,1,1);plot(x1); title('加入噪声后的信号波形'); subplot(2,1,2);
plot(f(1:600),abs(y1(1:600))); title('加入噪声后的信号频谱'); %高通 fp=3000; fs=2800; Fs=22050; rp=1;rs=100; wp=2*pi*fp/Fs; ws=2*pi*fs/Fs; Fs=Fs/Fs;
wap=2*tan(wp/2); was=2*tan(ws/2);
[N,wc]=cheb1ord(wp,ws,rp,rs,'s');[B,A]=cheby1(N,rp,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs); figure(4);
[h,w]=freqz(Bz,Az,512,Fs*22050); plot(w,abs(h));
title('切比雪夫高通滤波器'); xlabel('频率（HZ）'); ylabel('耗损（dB）'); grid on;
yd=filter(Bz,Az,x1); figure(5);
subplot(3,1,2); plot(yd);
title('滤波后的时域波形图'); ydd=fft(yd,800); subplot(3,1,3);
plot(f(1:600),abs(ydd(1:600)));
title('滤波后的频域波形图');
第 1页 共39页
语音信号的数字滤波处理(九)
ＩＶ切比雪夫带通滤波器
[x,Fs,bits]=wavread('C:\\Users\\acer\\Desktop\\ca9e36e190eb780c00b59c4459a6d5.wav'); x=x(:,1); figure(1); subplot(2,1,1); plot(x);
sound(x,Fs,bits);
title('语音信号时域波形图'); y=fft(x,3260);
f=(Fs/1630)*[1:1630]; subplot(2,1,2);
plot(f(1:1630),abs(y(1:1630))); title('语音信号频谱图'); %被污染的语音
t=0:length(x)-1;
zs=0.05*cos(2*pi*10000*t/22050); zs0=0.05*cos(2*pi*10000*t/22050000); figure(2); subplot(2,1,1) plot(zs0)
title('噪声信号波形'); zs1=fft(zs,1200);
sound(zs,Fs,bits); %回放噪音 subplot(2,1,2)
plot(f(1:600),abs(zs1(1:600))); title('噪声信号频谱'); %加噪声之后的语音 x1=x+zs';
%sound(x1,Fs,bits); %回放加入噪声后的语音 y1=fft(x1,1200); figure(3);
subplot(2,1,1);plot(x1); title('加入噪声后的信号波形'); subplot(2,1,2);
plot(f(1:600),abs(y1(1:600)));
title('加入噪声后的信号频谱');
fs1=500;fp1=800;fpu=1500;fsu=2000;Fs=22050' rp=1; rs=100;
wp=2*pi*[fp1,fpu]/Fs; ws=2*pi*[fs1,fsu]/Fs; Fs=Fs/Fs;
第 1页 共39页
语音信号的数字滤波处理(九)
wap=2*tan(wp/2); was=2*tan(ws/2);
[N,wc]=cheb1ord(wp,ws,rp,rs,'s'); [B,A]=cheby1(N,rp,wc,'s'); [Bz,Az]=bilinear(B,A,Fs);
[h,w]=freqz(Bz,Az,512,Fs*22050); figure(4); plot(w,abs(h));
title('切比雪夫带通滤波器'); xlabel('频率（HZ）'); ylabel('耗损（dB）'); grid on;
yd=filter(Bz,Az,x1); figure(5); subplot(2,1,1); plot(yd);
title('滤波后的时域波形图'); ydd=fft(yd,800); subplot(2,1,2);
plot(f(1:600),abs(ydd(1:600))); title('滤波后的频域波形图'); sound(yd,Fs,bits);
第 1页
共39页
语音信号的数字滤波处理(九) 附录Ｂ 加窗滤波器完整程序 I hamming窗低通滤波器完整程序
fs=22050;
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); y=y(1:32768)+0.5; figure(1) subplot(2,1,1); plot(abs(y)); xlabel('n'); ylabel('y');
title('原语音信号时域波形') y=y-0.5; Y=fft(y,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Y(1:16385))); xlabel('f'); ylabel('Y');
title('原语音信号频谱图') t=0:0.000045:0.90855;
x1=0.06*sin(2*pi*5500*t)+0.07*sin(2*pi*6000*t); x=[rot90(x1,2),rot90(x1,2)]; x=x(1:32768); z=x+y; figure(2) subplot(2,1,1); plot(z); xlabel('n'); ylabel('z');
title('被污染语音信号时域波形') Z=fft(z,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Z(1:16385))); xlabel('f'); ylabel('Z');
title('被污染语音信号频谱图')
wpz=2*1000*pi/fs;wsz=2*1200*pi/fs;rs=100; Bt=wsz-wpz;
N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2);
第 1页 共39页
语音信号的数字滤波处理(九)
wc=(wpz+wpz)/2/pi;
hn=fir1(N-1,wc,hamming(N)); figure(3)
freqz(hn,1,32678); figure(4)
hn=filter(hn,1,z); H=fft(hn,32768);
f=22050*(0:16382)/32768; sound(H,fs,nbits); subplot(2,1,2);
plot(f,abs(H(1:16383))); xlabel('f'); ylabel('H');
title('滤波后语音信号频谱图') subplot(2,1,1); hn=hn+0.5; plot(abs(hn)); xlabel('n'); ylabel('hn');
title('滤波后语音信号时域波形')
ＩＩhamming高通滤波器完整程序
fs=22050;
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); y=y(1:32768); figure(1) subplot(2,1,1); plot(y); xlabel('n'); ylabel('y');
title('原语音信号时域波形') Y=fft(y,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Y(1:16385))); xlabel('f'); ylabel('Y');
title('原语音信号频谱图') t=0:0.000045:0.90855;
x1=0.1*sin(2*pi*100*t)+0.2*sin(2*pi*220*t); x=[rot90(x1,2),rot90(x1,2)]; x=x(1:32768);
第 1页 共39页
语音信号的数字滤波处理(九)
z=x+y; figure(2) subplot(2,1,1); plot(z); xlabel('n'); ylabel('z');
title('被污染语音信号时域波形') Z=fft(z,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Z(1:16385))); xlabel('f'); ylabel('Z');
title('被污染语音信号频谱图')
wpz=2*3000*pi/fs;wsz=2*2800*pi/fs; rs=150; Bt=wpz-wsz;
N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2); wc=(wpz+wpz)/2/pi;
hn=fir1(N-1,wc,'high',hamming(N)); figure(3)
freqz(hn,1,32678); figure(4)
hn=filter(hn,1,z); H=fft(hn,32768);
f=22050*(0:16382)/32768; sound(H,fs,nbits); subplot(2,1,2);
plot(f,abs(H(1:16383))); xlabel('f'); ylabel('Yk');
title('滤波后语音信号频谱图') yk=ifft(hn,32768); subplot(2,1,1); plot(hn); xlabel('n'); ylabel('hn');
title('滤波后语音信号时域波形')
第 1页 共39页
语音信号的数字滤波处理(九)
ＩＶ hamming带通滤波器完整程序
fs=22050;
[y,fs,nbits]=wavread('C:\\Users\\acer\\Desktop\\邓紫棋-你把我灌醉.wav'); y=y(1:32768); figure(1) subplot(2,1,1); plot(y); xlabel('n'); ylabel('y');
title('原语音信号时域波形') Y=fft(y,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Y(1:16385))); xlabel('f'); ylabel('Y');
title('原语音信号频谱图') t=0:0.000045:0.90855;
x1=0.1*sin(2*pi*100*t)+0.2*sin(2*pi*220*t); x=[rot90(x1,2),rot90(x1,2)]; x=x(1:32768); z=x+y; figure(2) subplot(2,1,1); plot(z); xlabel('n'); ylabel('z');
title('被污染语音信号时域波形') Z=fft(z,32768);
f=22050*(0:16384)/32768; subplot(2,1,2);
plot(f,abs(Z(1:16385))); xlabel('f'); ylabel('Z');
title('被污染语音信号频谱图')
wpz=2*3000*pi/fs;wsz=2*2800*pi/fs; rs=150; Bt=wpz-wsz;
N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2); wc=(wpz+wpz)/2/pi;
第 1页 共39页
语音信号的数字滤波处理(九)
hn=fir1(N-1,wc,'high',hamming(N)); figure(3)
freqz(hn,1,32678); figure(4)
hn=filter(hn,1,z); H=fft(hn,32768);
f=22050*(0:16382)/32768; sound(H,fs,nbits); subplot(2,1,2);
plot(f,abs(H(1:16383))); xlabel('f'); ylabel('Yk');
title('滤波后语音信号频谱图') yk=ifft(hn,32768); subplot(2,1,1); plot(hn); xlabel('n'); ylabel('hn');
title('滤波后语音信号时域波形')
第 1页共39页
因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- huatuo9.cn 版权所有 赣ICP备2023008801号-1
违法及侵权请联系：TEL:199 18 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务