语音信号的数字滤波处理数字信号课程设计正文-大学论文.doc

一语音信号的数字滤波处理

目录

1 绪论……………………………………………………1
1.1 课题背景及目的………………………..……………………………1
1.2 设计方案及原理………………………..……………………………2
2 设计题第一部分………………………..………………………6
2.1 卷积的演示………………………..……………………6
2.2 采样定理的演示………………………..……………………12
3 设计题第二部分………………………..…………………15
3.1 设计要求………………………..…………………15
3.2 采集语音信号分析………………………..…………………16
4 设计数字滤波器………………………..…………………19
4.1 间接法设计IIR数字滤波器……………………..…………………19
4.1 直接法设计FIR数字滤波器……………………..…………………25
5 滤波分析……………………..…………………31
6 心得体会……………………..…………………35
7 参考文献……………………..…………………35
附录……………………..…………………36

1 绪论

1.1 设计目的
1.了解语音信号的产生、采集能绘制语音信号的频率响应曲线及频谱图 2. 学会用MATLAB对语音信号进行分析和处理
3. 学习设计滤波器的一般方法，通过课程设计树立正确的设计思想，提高分析问题、解决问题的能力
4.掌握用滤波器去除语音信号噪声的方法，观察去噪前后的语音信号。
5进行设计基本技能的训练，如查阅设计资料和手册、程序的设计、调试等。
1.2 设计要求
(1) 熟练掌握数字信号处理的基本概念基本理论。
(2) 熟悉离散信号和离散系统的时域特性。
(3) 掌握序列快速傅里叶变换（FFT）方法。
(4) 掌握滤波器的几种网络结构。
(5) 学会使用软件MATLAB，掌握其程序设计方法。
(6) 掌握利用MATLAB对语音信号进行频谱分析的方法。

1.2 课程设计方案
选择windows自带的“Windows XP 关机.wav”声音作为语音信号，产生噪声信号加到语音信号中，得到被污染的语音信号，经过频谱分析，分别设计巴特沃思滤波器以及Blackman窗滤波器（低通、带通、高通），对被污染的语音信号滤波，滤掉噪音信号得到符合要求的语音信号。
图2.1 总体方案设计流程图

图2.2 总体方案设计流程图（Blackman窗）

1.2.1 课程设计原理
设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H(z)具有指定的频率特性。数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应，（IIR）数字滤波器和有限长单位冲激响应（FIR）数字滤波器。本次数字信号处理课程设计我们小组选用巴特沃思数字滤波器（IIR）以及Blackman窗数字滤波器（FIR）。

1.2.2 IIR数字滤波器的结构与设计
IIR滤波器系统函数的极点可以在单位圆内的任何位置，实现IIR滤波器的阶次较低，所用的存储单元较少，效率高，又由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。设计IIR数字滤波器的方法主要有基于冲激响应不变法的IIR数字滤波器设计，基于双线性Z变换法的IIR数字滤波器设计，数字高通、带通及带阻IIR滤波器设计，基于MATLAB函数直接设计IIR数字滤波器。
1.2.3 基于双线性Z变换法的IIR数字滤波器设计
本次数字信号处理课程设计我们小组选用巴特沃思数字滤波器（IIR）以及Blackman窗数字滤波器（FIR）。由于 的频率映射关系是根据 推导的，所以使jΩ轴每隔2π/Ts便映射到单位圆上一周，利用冲激响应不变法设计数字滤波器时可能会导致上述的频域混叠现象。为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：
(1)s平面的整个jΩ轴仅映射为z平面单位圆上的一周。
(2)若G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的。
(3)这种映射是可逆的，既能由G(s)得到H(z)，也能由H(z)得到G(s)。
(4)如果G(j0)=1那么，双线性Z变换满足以上4个条件的映射关系，其变换公式为：
双线性Z变换的基本思路是，首先将整个s平面压缩到一条从–π/Ts到π/Ts的带宽为2π/Ts的横带里，然后通过标准的变换关系将横带变换成整个z平面上去，这样就得到s平面与z平面间的一一对应的单值关系。
在MATLAB中双线性Z变换可以通过bilinear函数实现，其调用格式为[Bz,Az]=bilinear(B,A,Fs)。
其中B,A为模拟滤波器传递函数G(s)的分子分母多项式的系数向量，而Bz,Az为数字滤波器的传递函数H(z)的分子分母多项式的系数向量。

1.2.4 数字高通、带通及带阻IIR滤波器设计
除了低通数字滤波器之外，实际中还常常需要高通、带通及带阻数字滤波器。这三种数字滤波器的设计步骤如下：
(1)将数字滤波器H(z)的技术指标ωp和ωs，通过Ω=tan(ω/2)转变为模拟滤波器G(s)的技术指标Ωp和Ωs，作归一化处理后，得到ηp=1，ηs=Ωs/Ωp；
(2)化解为模拟原型滤波器G(s)的技术指标；
(3)设计模拟原型滤波器G(p)；
(4)将G(p)转换为模拟滤波器的转移函数G(s)；
(5)将G(s)转换成数字滤波器的转移函数H(z)s=(z–1)(z+1)。
所谓原型滤波器是指归一化的低通滤波器。本节主要讨论通过IIR数字滤波器的原型转换设计法和IIR数字滤波器的直接设计方法来设计数字高通、带通及带阻滤波器其转换方法主要有3种一是直接由模拟低通滤波器转换成数字高通、带通或带阻滤波器二是先由模拟低通滤波器转换成模拟高通、带通或带阻滤波器然后再把它转换成相应的数字滤波器三是将模拟低通滤波器先转换成数字低通滤波器再通过变量代换变换成高通、带通或带阻滤波器。
1.3 设计Blackman窗数字滤波器的基本思想
1.3.1 FIR数字滤波器的结构与设计
FIR滤波器通常采用窗函数方法来设计。正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。
1.3.2 窗设计FIR数字滤波器的基本思想
窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后截取（加窗）它的脉冲响应得到线性相位和因果FIR滤波器。因此这种方法的重点是选择一个合适的窗函数和理想滤波器。我我们用表示理想的选频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽的低通滤波器由下式给定：

为了从得到一个FIR滤波器必须同时在两边截取。而要得到一个因果的线性相位滤波器，它的长度为N必须有：

这种操作叫做加窗，h(n)可以看作是与窗函数的乘积：

其中

根据的不同定义，可以得到不同的窗结构。
在频域中，因果FIR滤波器响应由和窗响应的周期卷积得到即：

2 设计题第一部分

2.1设计要求：
① 可输入任意两个序列x1(n)、x2(n)，指定x1(n)为自己的学号,例如x1(n)={2,0,1,2,5,7,0,5,0,1,1,4}。
x2(n)的内容和长度自选。例如x2(n)={ 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.0000}。
② 分别动态演示两个序列进行线性卷积x1(n)﹡x2(n)和圆周卷积x1(n)⊙x2 (n)的过程；要求分别动态演示翻转、移位、乘积、求和的过程；
③ 圆周卷积默认使用两个序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析；
④ 改变圆周卷积长度N，根据实验结果分析2类卷积的关系。
⑤ 在计算机操作系统下选一段声音文件(XP系统在“C:\WINDOWS\Media”)，读取文件取10ms的声音数据产生时域序列x1(n),序列内容自定义。利用x2(n)= {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}。利用FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比（验证时采用matlab子函数）。

2.2 卷积计算的演示
2.2.1 线性卷积的演示
演示程序：

第 43 页

n=-15:18;
M=17;
yn=zeros(1,19);
x1=[2,0,1,2,5,7,0,5,0,1,1,4];
xm1=[zeros(1,11),2,0,1,2,5,7,0,5,0,1,1,4,zeros(1,11)];
subplot(3,2,1);
stem(n,xm1)
xlabel('m')
ylabel('x1(m)')
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.0000];
subplot(3,2,2);
stem(xn2)
xlabel('m')
ylabel('x2(m)')
xm2=[fliplr(xn2) zeros(1,18)];
subplot(3,2,3);
stem(n,xm2)
xlabel('m')
ylabel('x2(N-m)')
title('n=0')
yn(1)=sum((xm1.*xm2)');
for N=1:17
xm3=[zeros(1,N) fliplr(xn2) zeros(1,M)];
M=M-1;
yn(N+1)=sum((xm1.*xm3));
end
xm3=[zeros(1,18) fliplr(xn2)]
yn(19)=sum((xm1.*xm3)');
subplot(3,2,4);
stem(yn)
xlabel('n')
ylabel('y(n)')
title('分步求卷积')
yn1=conv(x1,xn2);
subplot(3,2,5);
stem(yn1)
xlabel('n')
ylabel('y(n)')
title('调用函数直接求')

演示结果：

2.2.2圆周卷积程序如下：

clf;
x=[2,0,1,2,5,7,0,5,0,1,1,4];
y=[5,4,3,6,1,7,4,8,4,5,4,8,3,8,7,8];
subplot(2,1,1)
stem([x,zeros(1,28)]);
N1=length(x);
N2=length(y);
N=N1;
if N2>N1
N=N2;
end
M=input('请输入圆周卷积长度\n');
if M>=N
N=M;
end
if N==N1
x=x;
else
x=[x,zeros(1,N-N1)];
end
X=zeros(N,N);
x1=zeros(1,N);
pause(2)
for n=1:N
if n==1
x1(n)=x(n);
else
x1(n)=x(N-n+2);
end
end
X(1,1:N)=x1;
for i=2:N
pause(1)
for n=1:N
if n==1
X(i,n)=X(i-1,N);
else
X(i,n)=X(i-1,n-1);
end
end
subplot(2,1,1)
stem(X(i,1:N));
end
if N2<N
y=[y,zeros(1,N-N2)];
end
Y=y';
W=X*Y;
W=W';
n=1:N;
subplot(2,1,2)
stem(n,W);

当循环卷积区间长度N=40时，运行结果图：

当循环卷积区间长度N=20时，运行结果图：

分析总结：由以上分析可知，当循环卷积区间长度大于28时，不会出现混叠显现；当循环卷积长度小于28时，会出现混叠现象从而失真。

2.3利用FFT实现快速卷积
2.3.1 Matlab程序如下：
m1=wavread('d:\media\ding.wav',10);m1=m1';m1=m1(1,:);
m2=[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(421);stem(m1);ylabel('x1(n)');
subplot(422);stem(m2);ylabel('x2(n)');
y=conv(m1,m2);
subplot(423);stem(y);ylabel('y');
N=max(length(m1),length(m2));
X1=fft(m1,25);X2=fft(m2,25);
subplot(424);stem(X1);ylabel('X1');
subplot(425);stem(X2);ylabel('X2');
Y1=X1.*X2;
subplot(426);stem(Y1);ylabel('Y1');
Y2=fft(y);
subplot(427);stem(Y2);ylabel('Y2
波形图如下所示：

(2)编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应（n=14）：
采样程序如下所示：
t=0:0.0005:0.3;
xa=14*exp(-14*sqrt(2)*pi*t).*sin(14*sqrt(2)*pi*t);
subplot(6,1,1);plot(t,xa);title('origin xa')
t=0.015;
n=0:t:0.3;
pt=ones(0.3/t+1);
subplot(6,1,2);stem(n,pt);title('采样脉冲');
xaa=14*exp(-14*sqrt(2)*pi*n).*sin(14*sqrt(2)*pi*n);
subplot(6,1,3);stem(n,xaa);title('采样后信号');
y1=fft(xa,4096);
subplot(6,1,4);plot(abs(y1));title('采样频率2kHz');
y2=fft(xaa,4096);
subplot(6,1,5);plot(abs(y2));title('采样频率66Hz');
t=0.005;
n=0:t:0.3;
xa1=14*exp(-14*sqrt(2)*pi*n).*sin(14*sqrt(2)*pi*n);
y2=fft(xa1,4096);
subplot(6,1,6);plot(abs(y2));title('采样频率200Hz’）；

波形图如下图所示：

结果分析总结：当序列的长度为2^N时，可用FFT变换，再用两者频域相乘，逆变换成时域，计算速度比直接计算卷积快的多。当循环区间长度L>=N+M-1时，线性卷积与循环卷积相等，且无时域混叠现象。

3设计题第二部分（数字滤波器设计）

3.1设计要求：
① 利用Windows下的录音机或其他软件，进行语音信号的采集（*.wav）；
② 语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；
③ 产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；
④ 污染信号的频谱分析，画出被污染的语音信号时域波形和频谱；
⑤ 根据有关的频谱特性，采用间接法设计IIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个IIR滤波器）
a. 模拟滤波器类型：切比雪夫滤波器（低通、带通、高通）
b．总体要求：Matlab原程序+仿真波形+技术指标
⑥ 根据有关的频谱特性，采用直接法设计FIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个FIR滤波器）
a. 滤波器类型：Blackman窗（低通、带通、高通）
b．总体要求：Matlab原程序+仿真波形+技术指标+窗函数
⑦ 用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波；
⑧ 分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；
⑨ 回放语音信号。

3.2.1 采集语音信号，频谱分析
将电脑D:\media\sunshine.wav这一声音文件导入MATLAB软件中，并使用MATLAB绘出采样后的语音信号的时域波形和频谱图。
MATLAB程序如下：
clf;
[x,FS,bits]=wavread('d:\media\sunshine.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('语音信号频谱图');
运行程序得到语音信号时域和频谱图形如图所示:

3.2.2添加噪声，频谱分析：
MATLAB程序如下：
t=0:length(x)-1;
zs=0.1*cos(2*pi*10000*t/22050);
zs0=0.1*cos(2*pi*10000*t/22050000);
figure(2);
subplot(2,1,1); plot(zs0)
title('噪声信号波形');
zs1=fft(zs,1200);
subplot(2,1,2);
plot(f(1:600),abs(zs1(1:600)));
title('噪声信号频谱');
运行程序得到语音信号时域和频谱图形如图所示:

3.2.3添加噪声信号：
MATLAB程序如下：
y=x+zs';
figure(3);
subplot(2,1,1); plot(y)
title('被噪声污染的信号波形');
y1=fft(y,2000);
subplot(2,1,2);
plot(f(1:1000),abs(y1(1:1000)));
title('被噪声污染的信号频谱');
运行程序得到语音信号时域和频谱图形如图所示:

4设计数字滤波器

4.1采用间接法设计IIR数字滤波器
4.1.1 切比雪夫I型数字低通滤波器
先观察信号和噪声的频谱特性，选取合适的通带截止频率和阻带边界频率，故选取wp=2*pi*500;ws=2*pi*2000;

4.1.1.1 MATLAB程序如下：
wp=2*pi*500;ws=2*pi*1500;Rp=0.1;As=60;
[N1,wp1]=cheb1ord(wp,ws,Rp,As,'s');
[B1,A1]=cheby1(N1,Rp,wp1,'s');
fk=0:1500/512:1500;wk=2*pi*fk;
Hk=freqs(B1,A1,wk);
plot(fk/1000,20*log(abs(Hk)));
xlabel('频率（Hz）');ylabel('幅度（dB）')
axis([0,1.2,-160,5])
4.1.1.2 给被污染的信号加低通数字滤波器
wp=2*pi*500;ws=2*pi*2000;Rp=0.1;As=60;
[N1,wp1]=cheb1ord(wp,ws,Rp,As,'s');
[B,A]=cheby1(N1,Rp,wp1,'s');
fk=0:2000/512:2000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
figure(4);
plot(fk/1000,20*log(abs(Hk)));
xlabel('频率（Hz）');ylabel('幅度（dB）')
title('滤波器的频谱分析');
axis([0,1.2,-160,5]);
[Bz,Az]=impinvar(B,A,Fs); %脉冲响应不变法
M=filter(Bz,Az,y);
sound(M);%回放污染语音
figure(5);
subplot(2,1,1);plot(M);
title('滤波后被污染的信号波形');
M1=fft(M,2000);
subplot(2,1,2);
plot(f(1:1000),abs(M1(1:1000)));
title('滤波后被污染的信号频谱');

运行程序得到语音信号时域和频谱图形如图所示:

4.1.2设计数字带通滤波器
先将带通滤波器的指标装换为对应的低通数字滤波气的指标，设计相应的低通系统函数，对低通系统函数进行频率变换，得到带通数字滤波器

4.1.2.1 MATLAB程序：
wp=2*pi*[100,500];ws=2*pi*[10,800];Rp=1;As=20;
[N1,wp1]=cheb1ord(wp,ws,Rp,As,'s');
[B1,A1]=cheby1(N1,Rp,wp1,'s');
fk=0:1500/512:1500;wk=2*pi*fk;
Hk=freqs(B1,A1,wk);
plot(fk/1000,20*log(abs(Hk)));
xlabel('频率（Hz）');ylabel('幅度（dB）')
axis([0,1.2,-160,5])
4.1.2.2 给被污染的信号加带通数字滤波器
MATLAB程序：
wp=2*pi*[100,500];ws=2*pi*[10,800];Rp=1;As=20;
[N1,wp1]=cheb1ord(wp,ws,Rp,As,'s');
[B,A]=cheby1(N1,Rp,wp1,'s');
fk=0:2000/512:2000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
figure(4);
plot(fk/1000,20*log(abs(Hk)));
xlabel('频率（Hz）');ylabel('幅度（dB）')
title('带通滤波器的频谱分析');
axis([0,1.2,-160,5]);
[Bz,Az]=impinvar(B,A,Fs);
M=filter(Bz,Az,y);
sound(M); %回放污染语音
figure(5);
subplot(2,1,1);plot(M);
title('滤波后被污染的信号波形');
M1=fft(M,2000);
subplot(2,1,2);
plot(f(1:1000),abs(M1(1:1000)));
title('滤波后被污染的信号频谱');

4.1.2.3运行结果：

4.1.3 设计数字高通滤波器
4.1.3.1 设计指标：
先将高通滤波器的指标装换为对应的低通数字滤波气的指标，设计相应的低通系统函数，对低通系统函数进行频率变换，得到高通数字滤波器，本次取高通滤波器通带边界频率350Hz,阻带边界截止频率150Hz

4.1.3.2 MATLAB程序：

wp=2*pi*350;
ws=2*pi*150;
rp=3;rs=10;
fs=1000;
[n,wc]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(n,rs);
[a,b,c,d]=zp2ss(z,p,k);
[at1,bt1,ct1,dt1]=lp2hp(a,b,c,d,wc);
[at2,bt2,ct2,dt2]=bilinear(at1,bt1,ct1,dt1,fs);%双线性变换法
[num,den]=ss2tf(at2,bt2,ct2,dt2);
figure(4);
freqz(num,den,128,fs);
axis([150,400,-60,0])
grid on;
4.1.3.3 给被污染的信号加高通数字滤波器
wp=2*pi*350;
ws=2*pi*150;
rp=3;rs=10;
fs=1000;
[n,wc]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(n,rs);
[a,b,c,d]=zp2ss(z,p,k);
[at1,bt1,ct1,dt1]=lp2hp(a,b,c,d,wc);
[at2,bt2,ct2,dt2]=bilinear(at1,bt1,ct1,dt1,fs);%双线性变换法
[num,den]=ss2tf(at2,bt2,ct2,dt2);
figure(4);
freqz(num,den,128,fs);
axis([150,400,-60,0])
grid on;
M=filter(num,den,y);
sound(M,FS,bits); %回放污染语音
figure(5);
subplot(2,1,1);plot(M);
title('滤波后被污染的信号波形');
M1=fft(M,2000);
subplot(2,1,2);
plot(f(1:1000),abs(M1(1:1000)));
axis([0,1400,0,150])
title('滤波后被污染的信号频谱');

4.1.3.4 程序运行结果：

4.2采用直接法设计FIR数字滤波器
4.2.1设计Blackman窗-低通：
4.2.1.1 设计指标：
选取低通滤波器指标，通带边界截止频率500Hz，阻带边界截止频率1000Hz

4.2.1.2 Matlab程序：

fp=500;fs=1000;
wp=2*pi*fp/FS;
ws=2*pi*fs/FS;
Bt=ws-wp;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
[bz,az]=fir1(N-1,wc,blackman(N));
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([0,40,-0.2,0.5]);
title('h(n)波形');
subplot(2,1,2);
plot(w/pi*FS,20*log(abs(h)));
axis([0,3500,-400,2]);
title('低通-损耗函数曲线');
4.2.1.3 给被污染的信号加窗：
fp=500;fs=1000;
wp=2*pi*fp/FS;
ws=2*pi*fs/FS;
Bt=ws-wp;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
[bz,az]=fir1(N-1,wc,blackman(N));
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([0,40,-0.2,0.5]);
title('h(n)波形');
subplot(2,1,2);
plot(w/pi*FS,20*log(abs(h)));
axis([0,3500,-400,2]);
title('低通-损耗函数曲线');
M=filter(bz,az,y);
sound(M,FS,bits); %播放加窗后语音
figure(5);
subplot(2,1,1);plot(M);
title('滤波后被污染的信号波形');
M1=fft(M,2000);
subplot(2,1,2);
plot(f(1:500),abs(M1(1:500)));
title('滤波后被污染的信号频谱');

4.2.1.4 程序运行结果：

4.2.2设计Blackman窗-带通
4.2.2.1 设计指标：
选取带通指标参数：通带截止频率100Hz，500Hz，阻带截止频率30Hz,800Hz

4.2.2.2 Matlab程序：

wlp=100*2*pi/FS;wup=500*2*pi/FS;
wls=30*2*pi/FS;wus=1500*2*pi/FS;
Bt=wlp-wls;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=[(wls+wlp)/2/pi,(wus+wup)/2/pi];
[bz,az]=fir1(N-1,wc,blackman(N));
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([100,200,-0.2,0.5]);
title('h(n)波形');
subplot(2,1,2);
plot(w/pi*FS,20*log(abs(h)));
title('带通-损耗函数曲线');
4.2.2.3 给被污染的信号加窗：
wlp=100*2*pi/FS;wup=500*2*pi/FS;
wls=30*2*pi/FS;wus=1500*2*pi/FS;
Bt=wlp-wls;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=[(wls+wlp)/2/pi,(wus+wup)/2/pi];
[bz,az]=fir1(N-1,wc,blackman(N));
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([100,200,-0.2,0.5]);
title('h(n)波形');
subplot(2,1,2);
plot(w/pi*FS,20*log(abs(h)));
title('带通-损耗函数曲线');
%axis([0,1,-600,50]);
M=filter(bz,az,y);
sound(M,FS,bits); %播放加窗后语音
figure(5);
subplot(2,1,1);plot(M);
title('滤波后被污染的信号波形');
M1=fft(M,2000);
subplot(2,1,2);
plot(f(1:1000),abs(M1(1:1000)));
axis([0,2000,0,150]);
title('滤波后被污染的信号频谱');

4.2.2.4 程序运行结果：

4.2.3设计Blackman窗-高通：
4.2.3.1 设计指标：
选取带通指标参数：通带截止频率600Hz阻带截止频率400Hz

4.2.3.2 Matlab程序：
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([0,70,-0.2,1]);
title('h(n)波形');
subplot(2,1,2);
plot(w/pi*FS,20*log(abs(h)));
axis([300,1400,-300,100]);
title('高通-损耗函fp=600;fs=400;
wp=2*pi*fp/FS;
ws=2*pi*fs/FS;
Bt=wp-ws;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
[bz,az]=fir1(N-1,wc,'high',blackman(N));
数曲线');
4.2.3.3 给被污染的信号加窗：

fp=600;fs=400;
wp=2*pi*fp/FS;
ws=2*pi*fs/FS;
Bt=wp-ws;
N0=ceil(11*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
[bz,az]=fir1(N-1,wc,'high',blackman(N));
[h,w]=freqz(bz,az);
figure(4);
subplot(2,1,1);
plot([bz,az]);
axis([0,70,-0.2,1]);
title('h(n)波形');
subplot(2,1,2);