长安大学 数字信号处理课设.docx
《长安大学 数字信号处理课设.docx》由会员分享，可在线阅读，更多相关《长安大学 数字信号处理课设.docx（42页珍藏版）》请在冰豆网上搜索。
长安大学数字信号处理课设
《数字信号处理》课程设计报告
学院（部）信息工程学院
专业电子信息工程
班级24030901
学生姓名昆哥
学号24030901
目录
一、课程设计目的…………………………………………………3
二、课程设计原理…………………………………………………3
三、主要仪器及材料………………………………………………6
四、课程设计内容…………………………………………………6
1．语音信号采集……………………………………………6
2．语音信号频谱分析………………………………………6
3．设计数字滤波器和对信号滤波…………………………8
4．等波纹逼近法……………………………………………14
（1）低通滤波器…………………………………………14
（2）高通滤波器…………………………………………15
（3）带通滤波器…………………………………………17
（4）双线性变换法设计低通滤波器……………………19
（5）双线性变换法设计高通滤波器……………………24
（6）双线性变换法设计带通滤波器……………………30
五、课程设计结论…………………………………………………35
六、课程设计思考…………………………………………………36
七、体会与总结……………………………………………………36
八、参考文献………………………………………………………37
数字信号处理综合设计
、课程设计目的
1．学会MATLAB的使用，掌握MATLAB的程序设计方法；
2．掌握在Windows环境下语音信号采集的方法；
3．掌握数字信号处理的基本概念、基本理论和基本方法；
4．掌握MATLAB设计FIR和IIR数字滤波器的方法；
5．学会用MATLAB对信号进行分析和处理。
、课程设计原理
1．用窗函数法设计FIR滤波器
根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n），如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：
在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。
用窗函数wd（n）将hd（n）截断，并进行加权处理，得到
如果要求线性相位特性，则h（n）还必须满足：
根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。
要根据所设计的滤波特性正确选择其中一类。
例如，要设计线性相位低通特性可选择h（n）=h（N-1-n）一类，而不能选h（n）=-h（N-1-n）一类。
验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。
2.用双线性变换法设计IIR数字滤波器
脉冲响应不变法的主要缺点是产生频率响应的混叠失真。
这是因为从S平面到Ｚ平面是多值的映射关系所造成的。
为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。
也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。
这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1所示。
图1双线性变换的映射关系
为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现
（1）
式中,T仍是采样间隔。
当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。
将式
（1）写成
将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得
再将S1平面通过以下标准变换关系映射到Z平面
z=es1T
从而得到S平面和Z平面的单值映射关系为：
（2）
（3）
式
（2）与式（3）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换
式
（1）与式
（2）的双线性变换符合映射变换应满足的两点要求。
首先,把z=ejω，可得
（4）
即S平面的虚轴映射到Z平面的单位圆。
其次，将s=σ+jΩ代入式（4），得
因此
由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。
也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。
因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。
双线性变换法优缺点
双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。
这是因为S平面与Z平面是单值的一一对应关系。
S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。
这个关系如式（4）所示，重写如下：
上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图2所示。
由图2看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。
图2双线性变换法的频率变换关系
但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（4）及图2所示。
由于这种频率之间的非线性变换关系，就产生了新的问题。
首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3所示。
图3双线性变换法幅度和相位特性的非线性映射
对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。
也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。
、主要实验仪器及材料
微型计算机、Matlab6.5教学版、TC编程环境。
、课程设计内容
1．语音信号的采集
（1）利用windows下的录音机（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8.000KHz，8位，单声道）或其他软件，录制一段自己的话音，时间控制在1秒左右。
然后将音频文件保存为“wqk.wav”。
（2）在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。
2.语音信号的频谱分析
利用函数wavread对语音信号进行读取，将它赋值给某一变量。
首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，利用函数fft对信号进行快速傅里叶变换，得到信号的频谱特性：
程序如下：
z1=wavread（'F:
\wqk.wav'）;
plot（z1）;
图形分析如下：
z1=wavread（'F:
\wqk.wav'）;
y1=z1（1:
8192）;
Y1=fft（y1）;
n=0:
8191;
plot（n,Y1）;
3.设计数字滤波器和对信号滤波
（1）窗函数设计低通滤波器
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp=1000;fc=1200;As=100;Ap=1;Fs=8000;
wc=2*pi*fc/Fs;wp=2*pi*fp/Fs;
wdel=wc-wp;
beta=0.112*（As-8.7）;
N=ceil（（As-8）/2.285/wdel）;
wn=kaiser（N+1,beta）;
ws=（wp+wc）/2/pi;
b=fir1（N,ws,wn）;
figure
（1）;
freqz（b,1）;
x=fftfilt（b,z1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图形分析如下：
（2）窗函数设计高通滤波器
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp=2800;fc=3000;As=100;Ap=1;Fs=8000;
wc=2*pi*fc/Fs;wp=2*pi*fp/Fs;
wdel=wc-wp;
beta=0.112*（As-8.7）;
N=ceil（（As-8）/2.285/wdel）;
wn=kaiser（N,beta）;
ws=（wp+wc）/2/pi;
b=fir1（N-1,ws,'high',wn）;
figure
（1）;
freqz（b,1）;
x=fftfilt（b,z1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
（3）窗函数设计带通滤波器
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;Fs=8000;
wp1=2*pi*fp1/Fs;wc1=2*pi*fc1/Fs;wp2=2*pi*fp2/Fs;wc2=2*pi*fc2/Fs;
wdel=wp1-wc1;
beta=0.112*（As-8.7）;
N=ceil（（As-8）/2.285/wdel）;
ws=[（wp1+wc1）/2/pi,（wp2+wc2）/2/pi];
wn=kaiser（N+1,beta）;
b=fir1（N,ws,wn）;
figure
（1）;
freqz（b,1）
x=fftfilt（b,z1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
4、等波纹逼近法
（1）低通滤波器
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fs=8000;f=[1000,1200];
m=[1,0];
rp=1;rs=100;
dat1=（10^（rp/20）-1）/（10^（rp/20）+1）;dat2=10^（-rs/20）;
rip=[dat1,dat2];
[M,fo,mo,w]=remezord（f,m,rip,fs）;
hn=remez（M,fo,mo,w）;
figure
（1）;
freqz（hn）;
x=filter（hn,1,y1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,0.03]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;：
图如下：
（2）高通滤波器
程序：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fs=8000;f=[2800,3000];
m=[0,1];
rp=1;rs=100;
dat1=（10^（rp/20）-1）/（10^（rp/20）+1）;dat2=10^（-rs/20）;
rip=[dat2,dat1];
[M,fo,mo,w]=remezord（f,m,rip,fs）;
hn=remez（M,fo,mo,w）;
figure
（1）;
freqz（hn）;
x=filter（hn,1,y1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,0.03]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图形：
（3）带通滤波器
程序：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fs=8000;f=[1000,1200,3000,3200];
m=[0,1,0];
rp=1;rs=100;
dat1=（10^（rp/20）-1）/（10^（rp/20）+1）;dat2=10^（-rs/20）;
rip=[dat2,dat1,dat2];
[M,fo,mo,w]=remezord（f,m,rip,fs）;
hn=remez（M,fo,mo,w）;
figure
（1）;
freqz（hn）;
x=filter（hn,1,y1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,0.03]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图形：
（4）双线性变换法设计低通滤波器
选用butter
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp=1000;fc=1200;As=100;Ap=1;Fs=8000;
wc=2*fc/Fs;wp=2*fp/Fs;
[N,ws]=buttord（wc,wp,Ap,As）;
[b,a]=butter（N,ws）;
figure
（1）;
freqz（b,a,512,Fs）;
x=filter（b,a,z1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
选用cheby1
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp=1000;fc=1200;As=100;Ap=1;;Fs=8000;
wc=2*fc/Fs;wb=2*fp/Fs;
[n,wp]=cheb1ord（wc,wb,Ap,As）;
[b,a]=cheby1（n,Ap,wp）;
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,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
选用ellip
程序：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fp=1000;fc=1200;As=100;Ap=1;;Fs=8000;
wc=2*fc/Fs;wb=2*fp/Fs;
[n,wp]=ellipord（wc,wb,Ap,As）;
[b,a]=ellip（n,Ap,As,wp）;
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,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图形截图：
（5）双线性变换法设计高通滤波器
①选用butter
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fc=2800;fp=3000;As=100;Ap=1;Fs=8000;
wc=2*fc/Fs;wp=2*fp/Fs;
[N,ws]=buttord（wc,wp,Ap,As）;
[b,a]=butter（N,ws,'high'）;
figure
（1）;
freqz（b,a,512,Fs）;
x=filter（b,a,z1）;
X=fft（x,8192）;
figure
（2）;
subplot（2,2,1）;plot（abs（Y1））;axis（[0,4000,0,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
②选用cheby1
程序设计如下：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=fft（y1）;
fc=2800;fp=3000;As=100;Ap=1;Fs=8000;
wc=2*fc/Fs;wb=2*fp/Fs;
[n,wp]=cheb1ord（wc,wb,Ap,As）;
[b,a]=cheby1（n,Ap,wp,'high'）;
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,1.0]）;
title（'滤波前信号频谱'）;
subplot（2,2,2）;plot（abs（X））;axis（[0,4000,0,1.0]）;
title（'滤波后信号频谱'）;
subplot（2,2,3）;plot（z1）;
title（'滤波前信号波形'）;
subplot（2,2,4）;plot（x）;
title（'滤波后信号波形'）;
sound（x,fs,bits）;
图如下：
图18
选用ellip
程序：
clear;closeall
[z1,fs,bits]=wavread（'F:
\wqk.wav'）
y1=z1（1:
8192）;
Y1=ff
升级会员 