数字滤波器报告书(附MATLAB程序代码)

课程设计
课程名称__数字信号处理________ 题目名称__应用Matlab对语音信
__号进行频谱分析及滤波
学生学院__信息工程学院________ 专业班级__xx级通信工程xx班____ 学号__
学生姓名__ 指导教师__
2010年01月04日
题目名称应用Matlab对语音信号进行频谱分析及滤波
学生学院信息工程学院
专业班级xx级通信xx班
姓名xxxxxxxx
学号xxxxxxxx
一、课程设计目的
数字信号处理是一门以算法为核心，理论和实践性较强的学科。

是电子信息工程、通信工程专业、电子信息科学与技术专业的一门重要的专业技术基础课。

数字信号处理课程是在学习完数字信号处理的相关理论后，进行的综合性训练课程，其目的是：
1、使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现
方法；
2、增强学生应用Matlab语言编写数字信号处理的应用程序及分析、解决实际问
题的能力；
三、课程设计内容
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。

下面对各步骤加以具体说明。

2．1语音信号的采集
利用Windows下的录音机，录制一段自己的话音，时间在1 s内。

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

通过wavread函数的使用，学生很快理解了采样频率、采样位数等概念。

2．2语音信号的频谱分析
首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。

2．3设计数字滤波器和画出其频率响应
给出各滤波器的性能指标：
(1)低通滤波器性能指标fb＝1 000 Hz，fc＝1 200 Hz，As＝100 dB，Ap＝1 dB。

(2)高通滤波器性能指标fc＝4 800 Hz，fb＝5 000 Hz As＝100 dB，Ap＝1 dB。

(3)带通滤波器性能指标fb1＝1 200 Hz，fb2＝3 000 Hz，fc1＝1 000 Hz，fc2＝3 200 Hz，As＝100 dB，Ap＝1 dB。

用窗函数法和双线性变换法设计上面要求的3种滤波器。

在Matlab中，可以利用函数fir1设计FIR滤波器，可以利用函数butte，cheby1和ellip设计IIR滤波器；利用Matlab中的函数freqz画出各滤波器的频率响应。

2．4用滤波器对信号进行滤波
用自己设计的各滤波器分别对采集的信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。

2．5比较滤波前后语音信号的波形及频谱
要求学生在一个窗口同时画出滤波前后的波形及频谱。

2．6回放语音信号
在Matlab中，函数sound可以对声音进行回放。

其调用格式：sound(x，fs，bits)；可以感觉滤波前后的声音有变化。

2．7设计系统界面
设计处理系统的用户界面。

在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。

四、设计结果分析
双线性变换法分析：
双线性变换法的优点：不会有高于折叠频率的分量，避免了出现频率响应混叠响应。

双线性变换法的缺点：
1一个线性相位的模拟滤波器经双线性变换后就得到非线性相位的数字滤波器，不再保持原有的线性相位。

2、这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，几某一频率段的幅频响应近似等于某一常数，否则变换所产生的数字滤波器幅频响应相当于原模拟滤波器的幅频响应会有畸变。

窗函数法分析：
窗函数法的优点：有闭合形式的公式可循，设计实用简单。

窗函数法的缺点：通带、阻带的截止频率不易控制。

结果总结:只有选择的fc、fb与抽样频率Fs相匹配，才能得到与原声相近的声音输出。

五、设计心得体会
在做这个课程设计之前，我就已接触过MATLAB，学习过它的一些基本操作。

但这次课程设计中涉及到很多MATLAB的一些专门函数，对之前很少使用MATLAB 的我来说，是个不小的挑战。

在不久前的数字信号处理实验中，我又学习到一些用MATLAB设计数字滤波器的方法，知道了一些专门函数。

但即便如此，刚开始设计时还是觉得无处下手。

为此我又重新温习了一遍教材中有关IIR、FIR滤波器的内容，自己动手计算了一次实验手册上的IIR低通滤波器的有关参数。

在深感滤波器参数计算繁琐的情况下，我又上网搜索了一些用MATLAB设计数字滤波器的例程。

只是例程毕竟是别人写的，在没多少注解的情况下，要读懂并不容易。

为此，我再专门对一些感到疑惑的函数搜索一番，弄明白他们的使用形式。

最后，在众多例程中，我挑选了一个接近该课程设计要求的IIR低通滤波器例程作参考。

对其中的函数参数进行修改，得到了一个完全符合课程设计要求的IIR低通滤波器源程序。

再依葫芦画瓢，接连写出其它的滤波器源程序。

最终，我参考实验手册中例程框架，把已写出的几个滤波器源程序串在一起，得到了成品。

经过这次课程设计，我对IIR、FIR数字滤波器的设计步骤加深了认识。

同时，我也对MATLAB软件的强大感到吃惊。

在数学研究方面，MATLAB的确是首屈一指。

有MATLAB的强大函数库支持，原本由人手计算的诸多滤波器参数可交由电脑处理，节省了大量时间。

对MATLAB功能认识的进一步加深，我在今后的学习与工作中也能提高自己的效率。

六、参考资料
1、《数字信号处理及MATLAB实现》余成波杨如民等编著清华大学出版社出版
2、《数字信号处理教程》程佩青清华大学出版社出版
3.、《数字信号处理基础及MATLAB实现》周辉董正宏编著北京希望电子出版社出版
3、《MATLAB在数字信号处理中的应用》薛年喜主编清华大学出版社
七、源程序代码
clear;
close all;
[z1,fs,bits]=wavread('shiyan.wav')%读取存放于同一文件夹中的声音文件length=input('length=');%D是对原声取样长度
y1=z1(1:length);%y1是原声时域波形
Y1=fft(y1);%Y1是原声FFT变换
i=1;
K=menu('请选择滤波器类型','IIR低通滤波器','IIR高通滤波器',...
'IIR带通滤波器','FIR低通滤波器','FIR高通滤波器',...
'FIR带通滤波器','原声','IIR带阻滤波器', 'FIR带阻滤波器','Exit'); while(K~=10)
if K==1
%IIR低通滤波器
fb=input('fb=');
fc=input('fc=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=2*fc/Fs;
wb=2*fb/Fs;
[N,ws]=buttord(wb,wc,Ap,As)
[b,a]=butter(N,ws)
figure(i);i=i+1;
freqz(b,a,512,Fs);
title('IIR低通滤波器幅频相频图');
x=filter(b,a,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('IIR低通滤波后信号频谱');
subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('IIR低通滤波后信号波形'); sound(x,fs,bits);
elseif K==2
%IIR高通滤波器
fc=input('fc=');
fb=input('fb=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=2*fc/Fs;
wb=2*fb/Fs;
[N,ws]=buttord(wc,wb,Ap,As);
[b,a]=butter(N,ws,'high');
figure(i);i=i+1;
freqz(b,a,512,Fs);
title('IIR高通滤波器幅频相频图'); x=filter(b,a,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('IIR高通滤波后信号频谱'); subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('IIR高通滤波后信号波形'); sound(x,fs,bits);
elseif K==3
%IIR带通滤波器
fb1=input('fb1=');
fb2=input('fb2=');
fc1=input('fc1=');
fc2=input('fc2=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=[2*fc1/Fs,2*fc2/Fs];
wb=[2*fb1/Fs,2*fb2/Fs];
[N,ws]=buttord(wb,wc,Ap,As);
[b,a]=butter(N,ws,'bandpass');
figure(i);i=i+1;
freqz(b,a,512,Fs);
title('IIR带通滤波器幅频相频图'); x=filter(b,a,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('IIR带通滤波后信号频谱'); subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('IIR带通滤波后信号波形'); sound(x,fs,bits);
elseif K==4
%FIR低通滤波器
fb=input('fb=');
fc=input('fc=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=2*pi*fc/Fs;
wb=2*pi*fb/Fs;
wdel=wc-wb;%确定过渡带宽wdel
beta=0.1102*(As-8.7);%求凯泽窗参数beta
N=ceil((As-7.95)/2.286/wdel);%求凯泽窗阶数N wn=kaiser(N,beta);%凯泽窗函数
ws=(wb+wc)/2/pi;
B=fir1(N-1,ws,wn);
figure(i);i=i+1;
freqz(B,1);
title('FIR低通滤波器幅频相频特性');
x=fftfilt(b,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('FIR低通滤波后信号频谱');
subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('FIR低通滤波后信号波形');
sound(x,fs,bits);
elseif K==5
%FIR高通滤波器
fc=input('fc=');
fb=input('fb=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=2*pi*fc/Fs;
wb=2*pi*fb/Fs;
wdel=abs(wc-wb);%确定过渡带宽wdel
beta=0.1102*(As-8.7);%求凯泽窗参数beta
N=ceil((As-7.95)/2.286/wdel);%求凯泽窗阶数N wn=kaiser(N,beta)%凯泽窗函数
ws=(wb+wc)/2/pi;
B=fir1(N-1,ws,'high',wn);
figure(i);i=i+1;
freqz(B,1);
title('FIR高通滤波器幅频相频特性');
x=fftfilt(B,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('FIR高通滤波后信号频谱');
subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('FIR高通滤波后信号波形');
sound(x,fs,bits);
elseif K==6
%FIR带通滤波器
fb1=input('fb1=');
fb2=input('fb2=');
fc1=input('fc1=');
fc2=input('fc2=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wb1=2*pi*fb1/Fs;
wc1=2*pi*fc1/Fs;
wb2=2*pi*fb2/Fs;
wc2=2*pi*fc2/Fs;
wdel=wb1-wc1;%确定过渡带宽wdel
beta=0.1102*(As-8.7);%求凯泽窗参数beta
N=ceil((As-7.95)/2.286/wdel);%求凯泽窗阶数N wn= kaiser(N,beta);%凯泽窗函数
ws =[(wb1+wc1)/2/pi,(wb2+wc2)/2/pi];
B=fir1(N-1,ws,wn);
figure(i);i=i+1;
freqz(B,1);
title('FIR带通滤波器幅频相频特性');
x=fftfilt(B,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('FIR带通滤波后信号频谱')
subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('FIR带通滤波后信号波形');
sound(x,fs,bits);
elseif K==7
%原声
figure(i);
subplot(2,1,1);
plot(abs(Y1));
title('原始语音信号FFT频谱')
subplot(2,1,2);
plot(y1);
title('原始语音信号时域波形')
xlabel('时间 n');
ylabel('幅值');
sound(z1,fs,bits);
elseif K==8
%IIR带阻滤波器
fb1=input('fb1=');
fb2=input('fb2=');
fc1=input('fc1=');
fc2=input('fc2=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wc=[2*fc1/Fs,2*fc2/Fs];
wb=[2*fb1/Fs,2*fb2/Fs];
[N,ws]=buttord(wb,wc,Ap,As);
[b,a]=butter(N,ws,'stop');
figure(i);i=i+1;
freqz(b,a,512,Fs);
title('IIR带阻滤波器幅频相频图'); x=filter(b,a,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('IIR带阻滤波后信号频谱'); subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('IIR带阻滤波后信号波形'); sound(x,fs,bits);
elseif K==9
%FIR带阻滤波器
fb1=input('fb1=');
fb2=input('fb2=');
fc1=input('fc1=');
fc2=input('fc2=');
As=input('As=');
Ap=input('Ap=');
Fs=input('Fs=');
wb1=2*pi*fb1/Fs;
wc1=2*pi*fc1/Fs;
wb2=2*pi*fb2/Fs;
wc2=2*pi*fc2/Fs;
wdel=wb1-wc1;%确定过渡带宽wdel
beta=0.1102*(As-8.7);%求凯泽窗参数beta
N=ceil((As-7.95)/2.286/wdel);%求凯泽窗阶数N
wn= kaiser(N,beta);%凯泽窗函数
ws =[(wb1+wc1)/2/pi,(wb2+wc2)/2/pi];
B=fir1(N-1,ws,'stop',wn);
figure(i);i=i+1;
freqz(B,1);
title('FIR带阻滤波器幅频相频特性');
x=fftfilt(B,z1);
X=fft(x,length);
figure(i);i=i+1;
subplot(2,2,1);plot(abs(Y1));
title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));
title('FIR带阻滤波后信号频谱')
subplot(2,2,3);plot(z1);
title('滤波前信号波形');
subplot(2,2,4);plot(x);
title('FIR带阻滤波后信号波形');
sound(x,fs,bits);
end
K=menu('请选择滤波器类型','IIR低通滤波器','IIR高通滤波器',...
'IIR带通滤波器','FIR低通滤波器','FIR高通滤波器',...
'FIR带通滤波器','原声','IIR带阻滤波器', 'FIR带阻滤波器','Exit');
end
close all;。