中南大学现代信号处理课程设计2.docx
《中南大学现代信号处理课程设计2.docx》由会员分享，可在线阅读，更多相关《中南大学现代信号处理课程设计2.docx（26页珍藏版）》请在冰豆网上搜索。
中南大学现代信号处理课程设计2
一、课程设计目的：
1．全面复习课程所学理论知识，巩固所学知识重点和难点，将理论与实践很好地结合起来。
2.掌握信号分析与处理的基本方法与实现
3．提高综合运用所学知识独立分析和解决问题的能力；
4．熟练使用一种高级语言进行编程实现。
二、课程设计题目和具体设计（分题列出）；
（五）第五题
1.设计题目及设计要求
设有一信号
，设计各种IIR数字滤波器以实现：
1）低通滤波器，滤除
的成分，保留成分
2）高通滤波器，滤除
的成分，保留成分
3）带通滤波器，滤除
的成分，保留成分
4）带阻滤波器，滤除
的成分，保留成分
要求：
1）求出各个滤波器的阶数，设计各滤波器。
画出各滤波器的幅频和相频特性，计算滤波器的系统函数H（z）
2）画出滤波前后信号的时域、频域波形
2.设计思想及系统功能分析
本题主要考察四种基本IIR数字滤波器的设计，需要先对时域信号进行分析，确定各种滤波器的参数指标，然后根据matlab的一些固定调用格式对滤波器的阶数，截止频率等进行设计最后再生成最终滤波器，然后利用生成的滤波器对信号进行滤波。
3.关键部分的理论分析与计算
关键部分在于分析滤波器的参数指标时，要将频率都进行归一化，因为此处设计的是数字滤波器；另外，通带衰减一般小于3dB，阻带衰减一般大于30dB；设计带通、带阻滤波器时其截至频率分为上限截至频率和下限截至频率，书写方式为区间形式；另外，此处根据信号形式选择通、阻带截频时，应注意在所给信号频率的附近根据要求作适当的调整，不能刚好将所给的频率作为截止频率，因为滤波器不可能完全理想，为了保证信号能尽量无失真的滤波出来，就要通带范围稍微大一点。
4.程序源代码
%第五题
%
（1）设计数字低通滤波器
figure
（1）;
T=1
wp1=0.4;ws1=0.5;ap1=1;as1=30;%给出数字参数指标
[N1,wc1]=buttord（wp1,ws1,ap1,as1）;%计算阶数N和3Db处的频率
N1%在命令框中显示N
[B1,A1]=butter（N1,wc1）;%求巴特沃斯型数字滤波器
w=linspace（0,2*pi,1000）;%在0-2*pi内将w分成1000个点
h1=freqz（B1,A1,w）;%调用freqz函数，求解幅频响应
magh1=abs（h1）;%求出幅频特性
pha1=angle（h1）;%求出相频特性
subplot（321）;plot（w/pi,magh1,'r'）;
%画图
title（'低通滤波器幅频特性曲线'）;xlabel（'模拟频率'）;
subplot（322）;plot（w/pi,pha1,'r'）;
title（'低通滤波器相频特性曲线'）;xlabel（'模拟频率'）;
n=0:
72;%选取序列长度
w=pi*n/73;%分出与n相对应的w
xn=1+cos（pi*n/4）+cos（2*pi*n/3）;%原序列
subplot（323）;stem（n,xn,'.'）;%画出原序列
title（'xn的时域曲线'）;xlabel（'n'）;ylabel（'xn'）;
X=abs（xn*exp（-1*j*n'*w））;%原序列的序列傅里叶变换
subplot（3,2,4）;stem（w,X,'.'）;%画出序列傅里叶变换的曲线
title（'xn的频域曲线'）;xlabel（'w'）;ylabel（'X'）;
y1=filter（B1,A1,xn）;%利用filter函数恢复出滤波后的函数
subplot（325）;stem（n,y1,'.m'）;%画图
title（'低通滤波后的时域曲线'）;xlabel（'n'）;ylabel（'y1'）;
Y1=abs（y1*exp（-1*j*n'*w））;%恢复后的函数的序列傅里叶变换
subplot（3,2,6）;stem（w,Y1,'.m'）;%画图
title（'低通滤波后的频域曲线'）;xlabel（'w'）;ylabel（'Y1'）;
%
（2）设计数字高通滤波器
figure
（2）;
wp2=0.5;ws2=0.4;ap2=1;as2=30;
[N2,wc2]=buttord（wp2,ws2,ap2,as2）;
N2
[B2,A2]=butter（N2,wc2,'high'）;
w=linspace（0,2*pi,1000）;
h2=freqz（B2,A2,w）;
magh2=abs（h2）;
pha2=angle（h2）;
subplot（321）;plot（w/pi,magh2,'r'）;
title（'高通滤波器幅频特性曲线'）;xlabel（'模拟频率'）;
subplot（322）;plot（w/pi,pha2,'r'）;
title（'高通滤波器相频特性曲线'）;xlabel（'模拟频率'）;
n=0:
72;
w=pi*n/73;
xn=1+cos（pi*n/4）+cos（2*pi*n/3）;
subplot（323）;stem（n,xn,'.'）;
title（'xn的时域曲线'）;xlabel（'n'）;ylabel（'xn'）;
X=abs（xn*exp（-1*j*n'*w））;
subplot（3,2,4）;stem（w,X,'.'）;
title（'xn的频域曲线'）;xlabel（'w'）;ylabel（'X'）;
y2=filter（B2,A2,xn）;
subplot（325）;stem（n,y2,'.m'）;
title（'高通滤波后的时域曲线'）;xlabel（'n'）;ylabel（'y2'）;
Y2=abs（y2*exp（-1*j*n'*w））;
subplot（3,2,6）;stem（w,Y2,'.m'）;
title（'高通滤波后的频域曲线'）;xlabel（'w'）;ylabel（'Y2'）;
%（3）设计数字带通滤波器
figure（3）;
wp3=[0.1,0.4];ws3=[0.05,0.5];ap3=1;as3=30;
[N3,wc3]=buttord（wp3,ws3,ap3,as3）;
N3
[B3,A3]=butter（N3,wc3）;
w=linspace（0,2*pi,1000）;
h3=freqz（B3,A3,w）;
magh3=abs（h3）;
pha3=angle（h3）;
subplot（321）;plot（w/pi,magh3,'r'）;
title（'带通滤波器幅频特性曲线'）;xlabel（'模拟频率'）;
subplot（322）;plot（w/pi,pha3,'r'）;
title（'带通滤波器相频特性曲线'）;xlabel（'模拟频率'）;
n=0:
72;
w=pi*n/73;
xn=1+cos（pi*n/4）+cos（2*pi*n/3）;
subplot（323）;stem（n,xn,'.'）;
title（'xn的时域曲线'）;xlabel（'n'）;ylabel（'xn'）;
X=abs（xn*exp（-1*j*n'*w））;
subplot（3,2,4）;stem（w,X,'.'）;
title（'xn的频域曲线'）;xlabel（'w'）;ylabel（'X'）;
y3=filter（B3,A3,xn）;
subplot（325）;stem（n,y3,'.m'）;
title（'带通滤波后的时域曲线'）;xlabel（'n'）;ylabel（'y3'）;
Y3=abs（y3*exp（-1*j*n'*w））;
subplot（3,2,6）;stem（w,Y3,'.m'）;
title（'带通滤波后的频域曲线'）;xlabel（'w'）;ylabel（'Y3'）;
%（4）设计数字带阻滤波器
figure（4）;
wp4=[0.05,0.5];ws4=[0.1,0.4];ap4=1;as4=30;
[N4,wc4]=buttord（wp4,ws4,ap4,as4）;
N4
[B4,A4]=butter（N4,wc4,'stop'）;
w=linspace（0,2*pi,1000）;
h4=freqz（B4,A4,w）;
magh4=abs（h4）;
pha4=angle（h4）;
subplot（321）;plot（w/pi,magh4,'r'）;
title（'带阻滤波器幅频特性曲线'）;xlabel（'模拟频率'）;
subplot（322）;plot（w/pi,pha4,'r'）;
title（'带阻滤波器相频特性曲线'）;xlabel（'模拟频率'）;
n=0:
72;
w=pi*n/73;
xn=1+cos（pi*n/4）+cos（2*pi*n/3）;
subplot（323）;stem（n,xn,'.'）;
title（'xn的时域曲线'）;xlabel（'n'）;ylabel（'xn'）;
X=abs（xn*exp（-1*j*n'*w））;
subplot（3,2,4）;stem（w,X,'.'）;
title（'xn的频域曲线'）;xlabel（'w'）;ylabel（'X'）;
y4=filter（B4,A4,xn）;
subplot（325）;stem（n,y4,'.m'）;
title（'带阻滤波后的时域曲线'）;xlabel（'n'）;ylabel（'y4'）;
Y4=abs（y4*exp（-1*j*n'*w））;
subplot（3,2,6）;stem（w,Y4,'.m'）;
title（'带阻滤波后的频域曲线'）;xlabel（'w'）;ylabel（'Y4'）;
5.测试数据及必要的理论分析与比较
由figure1可以看出低通滤波器的幅频特性曲线和相频特性曲线，原信号的时域曲线和频域曲线以及滤波后的信号的时域曲线和频域曲线，通过对比，可以很清楚的看出该低通滤波器实现了题目所要求的滤波功能；
由figure2可以看出高通滤波器的幅频特性曲线和相频特性曲线，原信号的时域曲线和频域曲线以及滤波后的信号的时域曲线和频域曲线，通过对比，可以很清楚的看出该高通滤波器实现了题目所要求的滤波功能；
由figure3可以看出带通滤波器的幅频特性曲线和相频特性曲线，原信号的时域曲线和频域曲线以及滤波后的信号的时域曲线和频域曲线，通过对比，可以很清楚的看出该带通滤波器实现了题目所要求的滤波功能；
由figure4可以看出带阻滤波器的幅频特性曲线和相频特性曲线，原信号的时域曲线和频域曲线以及滤波后的信号的时域曲线和频域曲线，通过对比，可以很清楚的看出该带阻滤波器实现了题目所要求的滤波功能；
（六）第六题
1.设计题目及设计要求
（1）用Hanning窗设计一线性相位带通数字滤波器，要求：
N=15,
。
观察它的实际3dB和20dB带宽。
N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N变化的影响；
（2）分别改用矩形窗和Blackman窗，设计
（1）中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点；总结窗的不同长度和不同窗对滤波器的影响
2.设计思想及系统功能分析
本题主要考察利用窗函数法设计线性相位的数字滤波器，与之前的设计方法并无太大区别，先是确定滤波器的技术指标（归一化），然后通过调用窗函数的语句进行设计，只不过这里需要根据窗函数的不同选择相应的过渡带宽及选择信号的长度和滤波器的阶数；同样，带通滤波器的截止频率需要用区间形式写出；对于不同的窗函数，过渡带宽会不同。
3.关键部分的理论分析与计算
关键部分应该有确定截止频率（注意要归一化），还有选择不同的窗函数对应的过渡带宽；另外就是阶数和信号长度的关系：
信号长度=阶数+1，这是从资料书上看来的，所以每个滤波器的阶数都是按照这个等式来确定的。
4.程序源代码
%第六题
%
（1）用汉宁窗设计线性相位带通滤波器
%当长度N1为15时
N1=15;
w1=0.3*pi;w2=0.5*pi;
wc=[0.3,0.5];%归一化后的通带截频
B1=8*pi/N1;%过渡带宽
M1=N1-1;%单位脉冲响应的h（n）的长度N=M+1,此处的M1是阶数；
n1=0:
M1;
hn1=fir1（M1,wc,hanning（N1））;%调用firl函数求解单位脉冲响应
[h1,w1]=freqz（hn1,1,512）;%求其频率响应
magh1=abs（h1）;%求其幅频特性
pha1=angle（h1）;%求其相频特性
figure
（1）;
subplot（221）;plot（w1/pi,10*log（magh1））;%画出其幅频特性
grid;
title（'N=15时,汉宁窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（222）;plot（w1/pi,pha1）;%画出其幅频特性
title（'N=15时,汉宁窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
grid;
%当长度N2为45时
N2=45;
wc=[0.3,0.5];%归一化后的通带截频
B2=8*pi/N2;%过渡带宽
M2=N2-1;%单位脉冲响应的h（n）的长度N=M+1,此处的M1是阶数；
n2=0:
M2;
hn2=fir1（M2,wc,hanning（N2））;%调用firl函数求解单位脉冲响应
[h2,w2]=freqz（hn2,1,512）;%求其频率响应
magh2=abs（h2）;%求其幅频特性
pha2=angle（h2）;%求其相频特性
subplot（223）;plot（w2/pi,10*log（magh2）,'g'）;%画出其幅频特性
grid;
title（'N=45时,汉宁窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（224）;plot（w2/pi,pha2,'g'）;%画出其幅频特性
title（'N=45时,汉宁窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
grid;
%
（2）用矩形窗设计线性相位带通滤波器
%当长度N3为15时
figure
（2）;
N3=15;
wc=[0.3,0.5];
B3=4*pi/N3;%过渡带宽
M3=N3-1;
n3=0:
M3;
hn3=fir1（M3,wc,boxcar（N3））;
[h3,w3]=freqz（hn3,1,512）;
magh3=abs（h3）;
pha3=angle（h3）;
subplot（221）;plot（w3/pi,10*log（magh3））;
grid;
title（'N=15时,矩形窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（222）;plot（w3/pi,pha3）;
grid;
title（'N=15时,矩形窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
%当长度N2为45时
N4=45;
wc=[0.3,0.5];
B4=4*pi/N4;%过渡带宽
M4=N4-1;
n4=0:
M4;
hn4=fir1（M4,wc,boxcar（N4））;
[h4,w4]=freqz（hn4,1,512）;
magh4=abs（h4）;
pha4=angle（h4）;
subplot（223）;plot（w4/pi,10*log（magh4）,'g'）;
grid;
title（'N=45时,矩形窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（224）;plot（w4/pi,pha4,'g'）;
grid;
title（'N=45时,矩形窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
%（3）用Blackman窗设计线性相位带通滤波器
%当长度N3为15时
figure（3）;
N5=15;
wc=[0.3,0.5];
B5=12*pi/N5;%过渡带宽
M5=N5-1;
n5=0:
M5;
hn5=fir1（M5,wc,blackman（N5））;
[h5,w5]=freqz（hn5,1,512）;
magh5=abs（h5）;
pha5=angle（h5）;
subplot（221）;plot（w5/pi,10*log（magh5））;
grid;
title（'N=15时,Blackman窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（222）;plot（w5/pi,pha5）;
grid;
title（'N=15时,Blackman窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
%当长度N2为45时
N6=45;
wc=[0.3,0.5];
B6=12*pi/N6;%过渡带宽
M6=N6-1;
n6=0:
M6;
hn6=fir1（M6,wc,blackman（N6））;
[h6,w6]=freqz（hn6,1,512）;
magh6=abs（h6）;
pha6=angle（h6）;
subplot（223）;plot（w6/pi,10*log（magh6）,'g'）;
grid;
title（'N=45时,Blackman窗的幅频特性'）;xlabel（'w'）;ylabel（'magh'）;
subplot（224）;plot（w6/pi,pha6,'g'）;
grid;
title（'N=45时,Blackman窗的相频特性'）;xlabel（'w'）;ylabel（'pha'）;
5.测试数据及必要的理论分析与比较
由figure1可以看出，N=15时的汉宁窗设计的滤波器的通带比N=45时的通带宽，旁瓣的衰减也要小一些，所以，阶数越大，所设计出来的滤波器效果越好，通带稳定，衰减快。
同理，根据fugure2和figure3的矩形窗和Blackman窗也可以得出相同的结论。
通过对三个图不同窗函数设计的滤波器进行比较，可以看出：
矩形窗设计的滤波器通带最稳定，衰减最快；汉宁窗次之；Blackman窗滤波效果相对最差。
所以，矩形窗更能较好的实现滤波效果。
（七）选做题7
1.设计题目及设计要求
音乐信号处理：
1）获取一段音乐或语音信号，设计单回声滤波器，实现信号的单回声产生。
给出加入单回声前后的信号频谱。
2）设计多重回声滤波器，实现多重回声效果。
给出加入多重回声后的信号频谱。
3）设计均衡器，使得得不同频率的混合音频信号，通过一个均衡器后，增强或削减某些频率区域。
2.设计思想及系统功能分析
设计思想为：
先从网上下载一段格式为wav的音乐信号，选择一个较长的长度，设计一个滤波器，然后用该滤波器对音乐信号进行滤波，实质上就是将信号拉长后，与原信号就行一次叠加，便可设计出单回声滤波器；给滤波器加入一个阶次，让滤波器进行多次滤波然后再与原信号叠加，便可得到多重回声滤波器。
3.关键部分的理论分析与计算
设计单回声滤波器的关键在于该滤波器的系统函数的分子部分，要通过补零的方式将滤波器进行拉长，最好还要加入一个固定常数项；设计多重回声滤波器时，系统函数分子部分补的零根据阶数的不同会相应增加，最后的常数项也要变为相应的阶次幂，而且分母部分也应该进行补零。
4.程序源代码
%第七题
%
（1）单回声滤波器
%addpath（'D:
ProgameFiles/MATLAB/R2011a/bin'）;
[x,fs,bits]=wavread（'D:
\Prettyboy.wav',2^16）;%读出音频文件
%[x,fs,bits]=wavread（'D:
\yue.wav',2^16）;%读出音频文件
wavplay（x,fs）;%播放音频文件
fs
bits
pause
（1）;
a1=0.8;
R1=5000;
B1=[1,zeros（1,R1-1）,0.8];%系统函数分子
A1=[1];
y=filter（B1,A1,x）;%滤波器函数
wavplay（y,fs）;
pause
（1）;
x1=fft（x）;
figure
（1）;
subplot（321）;plot（abs（x1））;
title（'单回声前信号幅频响应'）;
subplot（322）;plot（angle（x1））;
title（'单回声前信号相频响应'）;
y1=fft（y）;%快速傅里叶变换
subplot（323）;plot（abs（y1）,'m'）;
title（'单回声后信号幅频响应'）;
subplot（324）;plot（angle（y1）,'m'）;
title（'单回声后信号相频响应'）;
[h1,w1]=freqz（B1,A1）;%求滤波器的幅频响应函数
subplot（325）;plot（abs（h1）,'g'）;
title（'单回声滤波器幅频响应'）;
y2=ifft（h1）;
subplot（326）;stem（abs（y2）,'.g'）;
title（'单回声滤波器冲激响应'）;
%多重回声滤波器
pause
（2）;
[x,fs,bits]=wavread（'D:
\Prettyboy.wav',2^16）;%读出音频文件
%[x,fs,bits]=wavread（'D:
\yue.wav',2^16）;%读出音频文件
wavplay（x,fs）;%播放音频文件
pause（0）;
N=5;
a2=0.8;
R2=8000;
B2=[1,zeros（1,N*R2-1）,-0.8^N];
A2=[1,zeros（1,R2-1）,-0.8];
z=filter（B2,A2,x）;
wavplay（z,fs）;
z1=fft（z）;
x1=fft（x）;
figure
（2）;
subplot（321）;plot（abs（x1））;
title（'多重回声前信号幅频响应'）;
subplot（322）;plot（angle（x1））;
title（'多重回声前信号相频响应'）;
subplot（323）;plot（abs（z1）,'m'）;
title（'多重回声后幅频响应'）;
subplot（324）;plot（angle（z1）,'m'）;
title（'多重回声后相频响应'）;
[h2,w2]=freqz（B2,A2）;
subplot（325）;plot（abs（h2）,'g'）;
title（'多重回声滤波器幅频响应'）;
z2=ifft（h2）;
subplot（326）;stem（abs（z2）,'.g'）;
title（'多重回声滤波器冲激响应'）;
5.测试数据及必要的理论分析与比较
由figure1可以看出，声音信号通过单回声滤波器后整个频率域上的幅值都发生了一些变化，变得比以前稍大一些，这是因为滤波后的信号是原信号与回声信号叠加起来的；同理，由figure2可以看出，声音信号通过多重回声滤波器后整个频率域上的幅值也都发生了一些变化，变得比以前更大。
通过对比figure1和figure2可以看出，单回声滤波器的幅频响应幅度比相应的多重回声滤波器的幅度小，因为多重回声滤波器是经过N次的单回声滤波器滤波后叠加的效果；而且单回声滤波器的冲激响应只有一个值，多重回声滤波器的
升级会员 