语音信号的数字滤波处理

 

摘要

本课程设计介绍了基于Matlab的对语音信号采集、处理及滤波器的设计，并使之实现的过程。运用课程中的基本概念、基本原理、基本分析方法，用Matlab进行数字语音信号处理，并阐述了课程设计的具体方法、步骤和内容。综合运用本课程的理论知识进行频谱分析以及滤波器设计，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。此外，系统实现了对数字信号处理的工作原理实施设计方法；用双线性变换法设计 IIR 数字滤波器，利用数字滤波器对信号进行滤波，对数字滤波器进行计算机仿真方法，并对设计结果加以分析。根据有关的频谱特性，直接法设计FIR数字滤波器，选择窗函数，设计带通、高通、低通滤波方法进行滤波。

 

关键词： IIR滤波器；FIR滤波器；切比雪夫滤波器；汉宁窗；

目录

1绪论………………………………………………………………………………………………………………2

2设计原理…………………………………………………………………………………………………………4

  2.1卷积运算……………………………………………………………………………………………………4

2.2采样定理………………………………………………………………………………………………………4

2.3用双线性变换法设计IIR数字滤波器……………………………………………… ………………………5

2.4用窗函数法设计FIR滤波器.……………………………………………………………………… ……… 6

3设计的前期准备..…………………………………………………………… …………………………… …8

3.1语音信号的采集…………………………………………… ………………………………… ……………8

3.2线性卷积……………………………………………………………………………… ………………… …9

3.3采样定理………………………………………………………………………………………………………12

4滤波器设计流程图………………………………………………………………………………………………14

5数字滤波器设计及其应用………………………………………………………………………………………15

5.1含噪声语音信号的频谱分析…………………………………………………………………………………15

5.2双线性变换法设计IIR滤波器.………………………………………………………………………………17

5.3窗函数设计FIR滤波器………………………………… ……………………………………………………25

6心得与体会………………………………………………………………………………………………………33

7参考文献…………………………………………………………………………………………………………34

1绪论

 

     数字滤波器可以在语音信号分析中对声音进行处理，可以滤出不要的噪声，使声音更加清楚。本设计通过对语音信号进行采集，对语音信号进行时域与频域的分析，然后给语音信号加上噪声，通过切比雪夫滤波器进行高通、低通、带通的滤波。通过汉宁窗对声音进行过滤。然后对声音进行回放，对比前后声音信号的差异。实现滤波功能。理论依据：根据设计要求分析系统功能，掌握设计中所需理论（采样频率、采样位数的概念，采样定理； 时域信号的FFT分析；数字滤波器设计原理和方法，各种不同类型滤波器的性能比较），阐明设计原理。信号采集：采集语音信号，并对其进行FFT频谱分析，画出信号的时域波形图和频谱图。构造受干扰信号并对其进行FFT频谱分析：对所采集的语音信号加入干扰噪声，对语音信号进行回放，感觉加噪前后声音的变化，分析原因，得出结论。并对其进行FFT频谱分析，比较加噪前后语音信号的波形及频谱，对所得结果进行分析，阐明原因，得出结论。 数字滤波器设计：根据待处理信号特点，设计合适数字滤波器，绘制所设计滤波器的幅频和相频特性。信号处理：用所设计的滤波器对含噪语音信号进行滤波。对滤波后的语音信号进行FFT频谱分析。画出处理过程中所得各种波形及频谱图。对语音信号进行回放，感觉滤波前后声音的变化。比较滤波前后语音信号的波形及频谱，对所得结果和滤波器性能进行频谱分析，阐明原因，得出结论。

2设计原理

2.1卷积运算

 

卷积和乘积运算在频域和时域是一一对应的，两个信号在时域的卷积可以转化为求两者在频域的乘积后再反变换，同理在频域的卷积等时域的乘积。而信号的频域求解有快速傅里叶FFT算法。

卷积与傅里叶变换有着密切的关系。利用一点性质，即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，能使傅里叶分析中许多问题的处理得到简化。

由卷积得到的函数f*g 一般要比f 和g 都光滑。特别当g 为具有紧支集的光滑函数，f 为局部可积时，它们的卷积f * g 也是光滑函数。利用这一性质，对于任意的可积函数f，都可以简单地构造出一列逼近于f 的光滑函数列，这种方法称为函数的光滑化或正则化。

卷积的概念还可以推广到数列、测度以及广义函数上去。

2.2采样定理

采样定理，又称香农采样定理，奈奎斯特采样定理，是信息论，特别是通讯与信号处理学科中的一个重要基本结论.E.T.Whittaker(1915年发表的统计理论),克劳德·香农与Harry Nyquist都对它作出了重要贡献。另外,V.A.Kotelnikov也对这个定理做了重要贡献。

　采样是将一个信号（即时间或空间上的连续函数）转换成一个数值序列（即时间或空间上的离散函数）。采样定理指出，如果信号是带限的，并且采样频率高于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。

带限信号变换的快慢受到它的最高频率分量的限制，也就是说它的离散时刻采样表现信号细节的能力是有限的。采样定理是指，如果信号带宽不到采样频率的一半（即奈奎斯特频率），那么此时这些离散的采样点能够完全表示原信号。高于或处于奈奎斯特频率的频率分量会导致混叠现象。大多数应用都要求避免混叠，混叠问题的严重程度与这些混叠频率分量的相对强度有关。

 

2.3用双线性变换法设计IIR数字滤波器

 

脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Ｚ平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=esT₁将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。

双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。

S平面上Ω与Z平面的ω成非线性的正切关系。

在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。

但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变。

对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。

 

2.4用窗函数法设计FIR滤波器

 

根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。

根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，h_d(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。

根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。

验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。

3设计前期准备

 

3.1语音信号的采集

1.语音信号的采集：利用windows下的c盘（C:\WINDOWS\Media \recycle，16位，双声道），，然后将音频文件保存“H:\recycle .wav”

2在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

3语音信号的频谱分析

首先画出语音信号的时域波形

z1=wavread(' H:\recycle .wav ');

plot(z1);

title('zhangwenbo   原始信号时域波形');

xlabel('数据采集单位 n');ylabel('amplitude');

图像输出如图1.1

图1.1语音信号时域波形

语音信号的频谱特性图1.2

z1=wavread('H:\recycle.wav');

plot(z1);

xlabel('Frequency（hz）');ylabel('Amplitude');

title('zhangwenbo  语音信号频谱特性');

 

图1.2语音信号频谱特性

 

3.2线性卷积

   任意两个序列x1(n)、x2(n)，指定x1(n)为自己的学号,例如x1(n)={2,0,0,7,8,4,2,5,0,1,2,3}。x2(n)的内容和长度自选进行线性卷积。

   线性系统的含义，就是，这个所谓的系统，带来的输出信号与输入信号的数学关系式之间是线性的运算关系。因此，实际上，都是要根据我们需要待处理的信号形式，来设计所谓的系统传递函数，那么这个系统的传递函数和输入信号，在数学上的形式就是所谓的卷积关系。卷积关系最重要的一种情况，就是在信号与线性系统或数字信号处理中的卷积定理。利用该定理，可以将时间域或空间域中的卷积运算等价为频率域的相乘运算，从而利用FFT等快速算法，实现有效的计算，节省运算代价。

线性卷积程序如下

clear

xls=[2,0,0,7,8,4,2,5,0,1,4,1];%输入X数组

lx=length(xls);

hls=[1,1,1,1];%输入h数组

lh=length(hls);

lmax=max(lx,lh);

if lx>lh    nx=0;nh=lx-lh;

elseif  lx<lh  nx=lh-lx;nh=0;

else  nx=0;nh=0;

end

dt=0.5;%运算时间间隔

lt=lmax;%取长者为补0长度为基准

x=[zeros(1,lt),xls,zeros(1,nx),zeros(1,lt)];%将X补锝与h同长，在两边补上同长度的0

tl=(-lt+1:2*lt)*dt;

h=[zeros(1,2*lt),hls,zeros(1,nh)];

hf=fliplr(h);%将h的左右翻转，记为hf

y=zeros(1,3*lt);%y数组初始化

for k=0:2*lt;

    p=[zeros(1,k),hf(1:end-k)];%hf向又循环移位

    y1=x.*p*dt;%输入和翻转移位脉冲函数相乘，在乘dt

    yk=sum(y1);%相加，等于积分

    y(k+lt+1)=yk;%将结果放入数组Y

    subplot(4,1,1);stem(tl,x)%用stem表示每次卷积结果

   

    axis([-lt*dt,2*lt*dt,min(x),max(x)]),hold on;

ylabel('x(t)');

     title('zhangwenbo');

    subplot(4,1,2);stem(tl,p)

    axis([-lt*dt,2*lt*dt,min(p),max(p)]),

    ylabel('h(k-t)');

    subplot(4,1,3);stem(tl,y1)

    axis([-lt*dt,2*lt*dt,min(y1),max(y1)+eps]),

    ylabel('s=x*h(k-t)');

    subplot(4,1,4);stem(k*dt,yk)

    axis([-lt*dt,2*lt*dt,floor(min(y)+eps),ceil(max(y)+eps)]),hold on;

    ylabel('y(k)=sum(s)*dt)');

    if k==round(0.8*lt)   disp('pause press any key to countinue'),pause

    else pause(1),

    end

end

线性卷积输出如图3.2

图3.2线性卷积 

3.3采样定理

采样位数即采样值或取样值，用来衡量声音波动变化的参数，采样频率是指在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。

采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精。

采样程序如下：

 

% 时域采样理论验证

Tp=64/1000;

%产生M长采样序列x(n)

% Fs=1000;T=1/Fs;

Fs=1000;T=1/Fs;

M=Tp*Fs;n=0:M-1;

A=41;alph=41*2^0.5*pi;omega=pi*41*2^0.5;

xnt=A*exp(-alph*n*T).*sin(omega*n*T);

Xk=T*fft(xnt,M);  %M点FFT[xnt)]

yn='xa(nT)';

subplot(3,2,1);tstem(xnt,yn);box on;title('(a) Fs=1000Hz');

k=0:M-1;fk=k/Tp;

subplot(3,2,2);plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=1000Hz');

xlabel('f(Hz)');ylabel('幅度');axis([0,180,0,1.2*max(abs(Xk))])

% Fs=300;T=1/Fs;

Fs=300;T=1/Fs;

M=ceil(Tp*Fs);n=0:M-1;

xnt=A*exp(-alph*n*T).*sin(omega*n*T);

Xk=T*fft(xnt,M);  %M点FFT[xnt)]

yn='xa(nT)';

subplot(3,2,3);tstem(xnt,yn);box on;title('(b) Fs=300Hz');

k=0:M-1;fk=k/Tp;

subplot(3,2,4);plot(fk,abs(Xk));title('(b) T*FT[xa(nT),Fs=300Hz');

xlabel('f(Hz)');ylabel('幅度');axis([0,180,0,1.2*max(abs(Xk))])

% Fs=200;T=1/Fs;

Fs=200;T=1/Fs;

M=ceil(Tp*Fs);n=0:M-1;

xnt=A*exp(-alph*n*T).*sin(omega*n*T);

Xk=T*fft(xnt,M);  %M点FFT[xnt)]

yn='xa(nT)';subplot(3,2,5);tstem(xnt,yn);box on;title('(c) Fs=200Hz');

k=0:M-1;fk=k/Tp;

subplot(3,2,6);plot(fk,abs(Xk));title('(c) T*FT[xa(nT),Fs=200Hz');

xlabel('f(Hz)');ylabel('幅度');axis([0,180,0,1.2*max(abs(Xk))])时域采样理论的验证程序运行如图3.3所示。由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。 

  图3.3时域采样

4滤波器设计流程图

 

本课程设计主要是对语音信号进行加噪声，滤波。首先要对声音进行采集，才可以进行信号处理。有了信号接下来要考虑的就是声音进行加噪声。可以插入一个白噪声或者脉冲信号。观察声音在频谱和时域在加噪声前后的变化，及加噪声与没加噪声的区别。可以通过对语音信号的时域波形，频谱进行比较。然后对语音信号进行频谱分析，对比。再通过MATLAB，根据有关的频谱特性，采用间接法设计IIR数字滤波器——切比雪夫滤波器进行低通、带通、高通滤波，看噪声过滤情况。还可以根据有关的频谱特性，采用直接法设计FIR数字滤波器——汉宁窗进行低通、带通、高通滤波，对比滤波情况。然后对滤波后波形进行频谱分析、最后在Matlab中,函数sound可以对声音进行回放。其调用格式：sound(x,fs,bits);可以感觉滤波前后的声音有变化。

录制或打开一个wav音频文件

开始

窗函数设计FIR滤波器

用双线性变换法设计IIR滤波器

频谱分析

声音回放对比

结束

加入噪声信号

对信号进行频谱时域分析

低通滤波

高通滤波

带通滤波

  

5数字滤波器设计及其应用

 

5.1污染信号的频谱分析

对信号进行采集后，往往有噪声。对于这些噪声应该怎么处理。我们可以用滤波器对信号进行滤波、对频谱进行分析。通过matlab的绘图功能进行比较，观察加噪声与不加噪声频谱的不同。

1、污染信号的频谱分析程序如下

Clear

[z1,fs,bits]=wavread('H:\recycle');

y1=z1(1:8192);

t=0:1/22050:8191*(1/22050);

z=sin(20000*pi*t);

m=y1+0.00000005*z';

u=fft(m,8192);

f=22050*(0:4096)/8192;

wp=5000/11025;ws=6000/11025;

[n,wn]=buttord(wp,ws,3,20);

[b,a]=butter(n,wn);

freqz(b,a,8192,22050);

s=filter(b,a,m);

subplot(2,2,2);

plot(f,abs(u(1:4097)));

title('张文波 滤波前信号频谱 ');

xlabel('Frequency（hz）');ylabel('Amplitude')

subplot(2,2,1);

plot(m(1:8192));

title('张文波 滤波前信号 '); xlabel('n');ylabel('Amplitude')

 

subplot(2,2,3);

plot(s(1:8192));

title('张文波 滤波后信号');

y=fft(s,8192);%对x进行256点的fft

subplot(2,2,4);

plot(f,abs(y(1:4097)));%画出频域内的信号

title('张文波 滤波后信号频谱 ');

 xlabel('Frequency（hz）');ylabel('Amplitude')

 

2、加噪声语音信号时域波形和频谱如图4.1

 

图5.1加噪声语音信号时域波形与频谱

 

5.2双线性变换法设计低通滤波器

      Chebyshev 滤波器则比Butterworth 滤波器的截至特性要好，在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。但通带处的幅值有振荡。对于数字滤波器而言，可以采用不同阶数逼近相应滤波器，滤波器性能还与滤波器的阶数有关，一般而言，阶数越高，则逼近越精确，但计算代价也随之上升，所以性能与代价总需要寻求一个平衡点。本设计用Chebyshev 滤波器。

4.2.1选用cheby1高通

在设计数字滤波器时，cheby1函数中的参数Wn是一个相对量，其定义区间

为Wn∈[0,1]，其中1对应于0.5fs，fs为采样频率（单位Hz）；在设计滤波器时，Wn采用真实频率，单位为Hz。 [b,a] = cheby1(n,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶[b,a] = cheby1(n,Rp,Wn,'ftype')用于设计数字高通和数字带阻滤波器，即‘ftype’ = high时，返回截止频率为Wn的数字高通滤波器；

 

程序设计如下：

 

%高通滤波器选用cheby1

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=22050;

wc=2*fc/Fs;wb=2*fp/Fs;% 设置指标参数

[n,wp]=cheb1ord(wc,wb,Ap,As);%计算切比雪夫1型模拟高通滤波器阶数和同

%带边界频率

[b,a]=cheby1(n,Ap,wp,'high');% 切比雪夫1型滤波器系统函数系数

%作图部分

figure(1);

freqz(b,a);

x=filter(b,a,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,4000,0,2.4]);

title('张文波  滤波前信号频谱');

xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,4000,0,1.2]);

title('滤波后信号频谱');

 xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')

sound(x,fs,bits);

 

图形分析如下图5.2.1.1 图5.2.1.2

 

 

                      图5.2.1.1

 

                      图5.2.1.2 cheby1高通滤波

 

5.2.2选用cheby1低通

cheby1函数中的参数Wn是一个相对量，其定义区间为Wn∈[0,1]，其中1对应于0.5fs，fs为采样频率（单位Hz）；在设计滤波器时，Wn采用真实频率，单位为Hz。[b,a] = cheby1(n,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶ChebyshevI型数字低通滤波器，通带内波纹为Rp。b、a分别为滤波器传递函数的分子和分母系数向量（降幂排列）。

程序设计如下：

%双线性变换法设计低通滤波器

%选用cheby1

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fp=1000;fc=1200;As=40;Ap=1; ;Fs=22050; % 设置指标参数

wc=2*fc/Fs;wb=2*fp/Fs; % 设置指标参数

[n,wp]=cheb1ord(wc,wb,Ap,As); %计算切比雪夫1型模拟低通滤波器阶数和同

%带边界频率

 [b,a]=cheby1(n,Ap,wp，’low’); % 切比雪夫1型滤波器系统函数系数

figure(1);

[h,w]=freqz(b,a);

plot(w,abs(h));

x=filter(b,a,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,2.0]);

title('张文波  滤波前信号频谱');

 xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,1000,0,2.0]);

title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后号波形'); xlabel('n');ylabel('Amplitude')

sound(x,fs,bits);

 

图形分析如下图5.2.2.1 图5.2.2.2

                            图5.2.2.1

图5.2.1.2 cheby1低通滤波

 

5.2.2选用cheby1带通

cheby1函数中的参数Wn采用真实频率，单位为Hz。[b,a] = cheby1(n,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶ChebyshevI型数字低通滤波器，当Wn为二元向量，即Wn=[W1 W2] (W1<W2)时，[b,a] = cheby1(n,Rp,Wn)返回一个2n阶数字带通滤波器，其通带为W1<ω< W2。

 

程序设计如下：

%带通滤波器选用cheby1

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;As=30;Ap=1; Fs=22050;

 wc=[2*fc1/Fs,2* fc2/Fs];wb=[2*fp1/Fs,2*fp2/Fs]; % 设置指标参数

 [n,wpo]=cheb1ord(wc,wb,Ap,As); %计算切比雪夫1型模拟带通滤波器阶数和同带边界频率

 [b,a]=cheby1(n,Ap,wpo); % 切比雪夫1型滤波器系统函数系数

%作图部分

figure(1);

freqz(b,a);

x=filter(b,a,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,3000,0,2.5]);

title('张文波——滤波前信号频谱');

xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,3000,0,2.5]);

title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')

sound(x,fs,bits);

 

图形分析如下图5.2.2.1 图5.2.2.2

图5.2.3.1

图5.2.3.2 cheby1带通滤波

 

5.3窗函数设计滤波器

汉宁窗又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是 3个 sine（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。这里用的就是汉宁窗。

5.3.1窗函数设计低通滤波器

选择一个适当的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。根据窗函数最小阻带衰减的特性，可采用汉宁窗提供大于30dB的衰减。本设计选择汉宁窗，其过渡带为8π/N，因此具有较小的阶次。

 

      程序设计如下：

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fp=1000 ;fc=1200  ;As=100 ;Ap=1 ;Fs=22050 ;

wp=2*fp/Fs;ws=2*fc/Fs;

DB=ws-wp;           %计算过渡带宽度

N0=ceil(6.2*pi/DB); %根据表7.2.2汉宁窗计算所需h(n)长度N0

N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数

wc=(wp+ws)/2/pi;    %计算理想高通滤波器通带截止频率(关于π归一化）

hn=fir1(N-1,wc,'low',hanning(N));  %调用fir1计算高通FIRDFh(n)

%以下是绘图部分

figure(1);

freqz(hn,1);

x=fftfilt(hn,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,300,0,2.0]);

title('张文波   滤波前信号频谱');

xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,300,0,0.5]);

title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')

 

图形分析如下图5.3.1.1 图5.3.1.2

 

图5.2.3.1

 

 

                            图5.2.3.2 汉宁窗低通滤波器

 

5.3.2窗函数设计高通滤波器

数字高通滤波器的设计分两种情况，即无相移和相移为±π/2两种情况。当无相移时，N为奇数时，所设计的FIR数字高通滤波器为Ⅰ型滤波器；当N为偶数时，为Ⅱ型滤波器。当相移为±π/2时，N为奇数时，所设计的FIR数字高通滤波器为Ⅲ型滤波器。根据窗函数最小阻带衰减的特性，选择汉宁窗可达到 44dB的最小阻带衰减，其过渡带为6.2π/N，因此具有较小的阶次。

 

程序设计如下：

%窗函数设计高通滤波器

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=22050;

wp=2*fp/Fs;ws=2*fc/Fs;

DB=wp-ws;           %计算过渡带宽度

N0=ceil(6.2*pi/DB); %根据表7.2.2汉宁窗计算所需h(n)长度N0

N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数

wc=(wp+ws)/2/pi;    %计算理想高通滤波器通带截止频率(关于π归一化）

hn=fir1(N-1,wc,'high',hanning(N));  %调用fir1计算高通FIRDFh(n)

%以下是绘图部分

figure(1);

freqz(hn,1);

x=fftfilt(hn,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,800,0,1.0]);

title('张文波  滤波前信号频谱');

xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,800,0,1.0]);

title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')

sound(x,fs,bits);

 

图形分析如下图5

 图5.3.2汉宁窗低通滤波器

 

5.3.3窗函数设计带通滤波器

 

数字带通滤波器的设计也分两种情况，即无相移和相移为±π/2两种情况。当无相移时，N为奇数时，所设计的FIR数字带通滤波器为Ⅰ型滤波器；当N为偶数时，为Ⅱ型滤波器。当相移为±π/2时，N为奇数时，所设计的FIR数字带通滤波器为Ⅲ型滤波器；当N为偶数时，为Ⅳ型滤波器。汉宁窗通带纹波和阻带纹波均满足设计要求。这里选汉宁窗。

程序设计如下：

%窗函数设计带通滤波器

clear;close all

[z1,fs,bits]=wavread('H:\recycle.wav')

y1=z1(1:8192);

Y1=fft(y1);

fp1=1200 ;fp2=3000 ;fc1=1000  ;fc2=3200 ;As=100 ;Ap=1 ;Fs=22050 ;

wp1=2*pi*fp1/Fs; wc1=2*pi*fc1/Fs; wp2=2*pi*fp2/Fs; wc2=2*pi*fc2/Fs;

wdel=wp1-wc1;       %计算过渡带宽度

N0=ceil(6.2*pi/wdel); %根据表7.2.2汉宁窗计算所需h(n)长度N0

N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数

wc=[(wp1+wc1)/2/pi,(wp2+wc2)/2/pi];    %计算理想带通滤波器通带截止频率(关于π归一化）

hn=fir1(N-1,wc,hanning(N));  %调用fir1计算高通FIRDFh(n)

%以下是绘图部分

figure(1);

freqz(hn,1);

x=fftfilt(hn,z1);

X=fft(x,8192);

figure(2);

subplot(2,2,1);plot(abs(Y1));axis([0,2000,0,3.0]);

title('张文波   滤波前信号频谱');

xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,2);plot(abs(X));axis([0,2000,0,3.0]);

title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');

subplot(2,2,3);plot(z1);

title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')

subplot(2,2,4);plot(x);

title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')

sound(x,fs,bits);

 

图形分析如下图5.3.3.1图5.3

 图5.3.2 

图5.3.2汉宁窗带通滤波器

 

 

6心得与体会

这次课程设计是我学习最多的一次设计，在贺老师的指导下，我们每天都有自己的任务，对我们的学习更有针对性，当我们遇到问题时，在贺老师的指导下可以很好的解决，这是以前的设计中从来没有的。了我对对内容的理解、吸收有点困难。但是通过一周的学习，让我重新看了数字信号处理，matlab等方面的书籍，更加上对列子实验的重新做过，是我进步很大。再者老师给我们每个同学布置不同的课题，使我们坚信这次课设只能看自己。在自己的学习、老师的帮助下。终于将课程设计做完。

    其实这个设计只要我们掌握了双线性变换法分析：双线性变换的主要优点、双线性变换法的缺点。窗函数法分析相位响应有严格的线性,不存在稳定性问题, 设计简单这些规律就简单了。另外一定要熟悉matlab这个软件。

课设的几天是很充实的，脑子在不停的运转，感觉学习理论一定要于实践相结合才能学到东西，不要只停留在书本。这几天学到很多东西.大家一起发现问题,讨论问题,解决问题.留给自己的就是所需要的知识,这种很互动的学习方式让我们实在是受益匪浅啊!

虽然得对着电脑半天,一天的做程序，有点累也有些晕，可是当看到劳动成果时，真是别有一番滋味在心头啊！世上无难事，只怕有心人，的确如此。做完这个程序最大的收获就是只要你肯做就一定可以成功。

做完这个课程设计，我们提高了自信心，我突然发现程序其实不是那么难，多学习重视好的，因为享受劳动成果的滋味实在很美妙啊！忙碌了几天，在大家的共同努力下，我们总算将此数字滤波器设计出来。尽管用了很长时间，但仍然很高兴，因为在设计的过程中，让我了解到要设计一个型程序，查找资料是至关重要的，这过程艰辛，但只要你持之以恒，成功指日可待。

通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才是真正的知识，才能提高自己的实际动手能力和独立思考的能力。

 

 

 

7参考文献

 

[1]丁玉美等. 数字信号处理（第二版[M].西安电子科技大学出版社.

[2]王创新,文卉.数字信号处理试验指导书[M].长沙理工大学印刷（内部使用）.

[3]陈怀琛. 数字信号处理及其MATLAB实现[M].电子工业出版社.

[4]陈怀琛等. MATLAB及在电子信息课程中的应用[M].电子工业出版社.

[5] A.V.奥本海姆，R.W.谢弗.数字信号处理[M].北京：科学出版社.

[6]胡广书. 数字信号处理——理论、算法与实现（第二版）[M]北京：电子工业出版社.