数字信号处理课程设计.docx
《数字信号处理课程设计.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计.docx（33页珍藏版）》请在冰点文库上搜索。
数字信号处理课程设计
《数字信号处理》课程设计
设计题目：
基于MATLAB的音乐信号处理和分析
一、课程设计的目的
本课程设计通过对音乐信号的采样、抽取、调制、解调等多种处理过程的理论分析和MATLAB实现，使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；使学生掌握的基本理论和分析方法知识得到进一步扩展；使学生能有效地将理论和实际紧密结合；增强学生软件编程实现能力和解决实际问题的能力。
二、课程设计基本要求
1学会MATLAB的使用，掌握MATLAB的基本编程语句。
2掌握在Windows环境下音乐信号采集的方法。
3掌握数字信号处理的基本概念、基本理论和基本方法。
4掌握MATLAB设计FIR和IIR数字滤波器的方法。
5掌握使用MATLAB处理数字信号、进行频谱分析、设计数字滤波器的编程方法。
三、课程设计内容
1、音乐信号的音谱和频谱观察
使用windows下的录音机录制一段音乐信号或采用其它软件截取一段音乐信号（要求：
时间不超过5s、文件格式为wav文件）
①使用wavread语句读取音乐信号，获取抽样率；（注意：
读取的信号是双声道信号，即为双列向量，需要分列处理）；
②输出音乐信号的波形和频谱，观察现象；
③使用sound语句播放音乐信号，注意不同抽样率下的音调变化，解释现象。
查找help：
Wavread格式说明：
[y,fs,b]=wavread（‘语音信号’），采样值放在向量y中，fs表示采样频率（hz），b表示采样位数。
【matlab程序如下】
clc
closeall
clearall
[y,fs]=wavread（'一生有你'）;%读取歌曲
size（y）
y1=y（:
1）;%1声道
y2=y（:
2）;%2声道
N=length（y）;%长度
n=0:
N-1;
t=n/fs;%t=nT
w=2*n/N;%2pi在长度N上的平均是个序列
f1=fft（y1）;
f2=fft（y2）;%傅里叶变换
figure%画图像
subplot（2,2,1）,plot（t（1:
1000）,y1（1:
1000））
title（'1声道时域图'）
xlabel（'t'）
ylabel（'y1'）
subplot（2,2,3）,plot（w,abs（f1）*2/N）
title（'1声道频域图'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（2,2,2）,plot（t（1:
1000）,y2（1:
1000））
title（'2声道时域图'）
xlabel（'t'）
ylabel（'y2'）
subplot（2,2,4）,plot（w,abs（f2）*2/N）
title（'2声道频域图'）
xlabel（'w'）
ylabel（'Y2（w）'）
wavplay（y1,fs）
wavplay（y1,fs/2）慢放wavplay（y1,fs*2）快放
wavylay（y2,fs*2）快放wavplay（y2,fs/2）慢放
【程序运行结果如下图】：
分析：
通过观察音乐信号的波形和频谱可知所选取的音乐信号频谱集中在0~0.5pi之间，同时抽样频率为fs=44000
2、音乐信号的抽取（减抽样）
①观察音乐信号频率上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；
②输出减抽样音乐信号的波形和频谱，观察现象，给出理论解释；
③播放减抽样音乐信号，注意抽样率的改变，比较不同抽取间隔下的声音，解释现象。
理论基础：
时域抽样定理：
一个频谱受限的信号f（t），如果频谱只占据-wm~+wm的围，
则信号f（t）可以用等间隔的抽样值唯一的表示。
而抽样间隔必须不大于1/（2*fm）.
频域抽样定理：
一个频谱受限的信号f（t），它集中在-tm~+tm的时间范围内，
若在频域中以不大于1/（2*tm）的频率间隔对f（t）的频谱F（w）进行抽样，则抽样后
频谱F1（w）可以唯一的表示原信号。
【Matlab程序如下】：
%原信号的频率上限为0.5pi
clearall;closeall
[y,fs,bits]=wavread（'一生有你'）;
y1=y（:
1）;%取一频道信号
f1=fft（y1）;
N=length（y1）;
%减抽样$$$$$$$$$$$$$$$$$$$减抽样使抽样点数减少，会使栅栏效应更严重
D1=2;D2=16;
n=0:
N-1;
t=n/fs;%t=nT
yd1=y1（1:
D1:
N）;%2倍减抽样
fyd1=fft（yd1,N）;
yd2=y1（1:
D2:
N）;
fyd2=fft（yd2,N）;%16倍减抽样
w=2*n/N;%2pi在长度N上的平均
figure
subplot（2,3,1）;plot（t（1:
1000）,y1（1:
1000））;
title（'原信号时域图'）
xlabel（'t'）
ylabel（'y1'）
subplot（2,3,4）;plot（w,abs（f1））;
title（'原信号频域图'）
xlabel（'t'）
ylabel（'y1'）
subplot（2,3,2）;plot（t（1:
1000）,yd1（1:
1000））;
title（'2倍减抽样后的时域图'）
xlabel（'t'）
ylabel（'yd1'）
subplot（2,3,5）;plot（w,abs（fyd1））;
title（'2倍减抽样后的频域图'）%出现栅栏效应
xlabel（'t'）
ylabel（'fyd1'）
subplot（2,3,3）;plot（t（1:
1000）,yd2（1:
1000））;
title（'10倍减抽样后的时域图'）
xlabel（'t'）
ylabel（'yd2'）
subplot（2,3,6）;plot（w,abs（fyd2））;
title（'16倍减抽样后的频域图'）%栅栏效应更明显
xlabel（'t'）
ylabel（'fyd2'）
wavplay（yd1,fs）;
wavplay（yd1,fs/D1）;
wavplay（yd2,fs）;
wavplay（yd2,fs/D2）;
%D倍减抽样后相当于原抽样频率缩小了D倍，故仍按原抽样频率播放相当于加快播放
【程序运行结果如下图】：
分析：
通过观察两种不同抽样间隔（2倍频和16倍频）下的音乐信号可知，当采用较大的抽样间隔对音乐信号进行抽样时，频谱发生了混叠，而采用较小的抽样间隔对音乐信号进行抽样时，频谱并未发生混叠。
这是因为，抽样时频谱发生混叠的条件是fs<2fh，即抽样频率小于信号频谱的最高频率。
当采用较大的抽样间隔时抽样频率时fs<2fh，所以发生混叠，而采用较小的抽样间隔时抽样频率时fs>=2fh，则不会发生混叠。
当我们播放不同抽样间隔下的音乐信号时，会发现大抽样的音乐信号会伴有杂音并且声音低沉，而小抽样的音乐信号和原有的音乐信号几乎无差别，这间接证明了我们以上理论分析的正确性。
3、音乐信号的AM调制
①观察音乐信号频率上限，选择适当调制频率对信号进行调制（给出高、低两种调制频率）；
②输出调制信号的波形和频谱，观察现象，给出理论解释；
③播放调制音乐信号，注意不同调制频率下的声音，解释现象。
理论基础：
信号的调制过程就是将信号频谱搬移到任何所需的较高频率范围。
调制的实质是把各种信号的频谱搬移，使它们互不重叠的占据不同的频率范围，也即信号分别托付于不同频率的载波上。
具体的调制原理推导在此不再叙述，仅将结论列出：
f（t）=g（t）*cos（w0*t）。
由此将信号g（t）的频谱搬移到（2*n+1）w附近，同时音乐信号的频谱幅度变为原来的1/2。
如果信号的最高频谱wh超过了ws/2，则各周期延括分量产生频谱的交叠，称为是频谱的混叠现象。
根据奈奎斯特定律可知，若希望频谱不会发生混叠，则fs>=fh。
【Matlab程序如下】：
%原信号的频率上限为0.5pi
clc
clearall
[y,fs]=wavread（'一生有你'）;
y1=y（:
1）;
N=length（y）;
n=0:
（N-1）;
%低调制频率
y2=cos（n*0.5*pi）;%调制频率为0.5*pi；
y3=y1.*y2';%调制后的信号；
f1=fft（y1）;%对原信号做fft变换；
f3=fft（y3）;%对调制后的信号做fft变换；
t=n/fs;%t=nT
N=length（y）;
w=2*n/length（y）;
%高调制频率
y4=cos（n*0.8*pi）;
y5=y1.*y4';
f5=fft（y5）;
%图形显示
figure
subplot（2,3,1）,plot（t,y1）
title（'原信号时域图'）
xlabel（'t'）
ylabel（'y1'）
subplot（2,3,4）,plot（w,abs（f1）*2/N）
title（'原信号频域图'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（2,3,2）,plot（t,y3）
title（'低调制后的时域图'）
xlabel（'t'）
ylabel（'y3'）
subplot（2,3,5）,plot（w,abs（f3）*2/N）
title（'低调制后的频域图'）
xlabel（'w'）
ylabel（'Y3（w）'）
subplot（2,3,3）,plot（t,y5）
title（'高调制后的时域图'）
xlabel（'t'）
ylabel（'y5'）
subplot（2,3,6）,plot（w,abs（f5）*2/N）
title（'高调制后的频域图'）
xlabel（'w'）
ylabel（'Y5（w）'）
%播放声音
wavplay（y1,fs）%播放调制前的原信号；
wavplay（y3,fs）%播放低调制后的信号；
wavplay（y5,fs）%播放高调制后的信号
【程序运行结果如下图】：
分析：
通过观察原音乐信号的频谱图可知音乐信号的频率上限是0.5pi，但是为了方便以后的计算，在此将0.3pi后的信号舍去，即默认音乐信号的频率上限是0.3pi。
取高频调制频率0.5pi，低频调制频率0.3pi对音乐信号进行调制。
通过观察不同调制频率下的频谱图可以发现高频调制的音乐信号频谱发生了混叠，而采用低频调制的音乐信号频谱并未发生混叠。
这是因为当采用高频调制（0.5pi）时，频谱被搬移到（2*n+1）*0,5pi，n=0.1.2.3…..附近，此时高频调制频率高于原信号的频率上限，故发生了频谱混叠。
同理，当采用低频调制（0.3pi）时，未发生频谱混叠。
播放不同调制频率下的音乐信号，可以发现当采用低频调制时，音乐信号比原信号的声音低了很多，但是没有杂音；采用高频调制时，音乐信号比原信号的声音低了很多的同时还伴随有杂音。
这是因为低频调制没有发生混叠，调制后的音乐信号频谱幅度为原音乐信号的1/2，而高频调制发生了混叠。
4、AM调制音乐信号的同步解调
①设计巴特沃斯IIR滤波器完成同步解调；观察滤波器频率响应曲线；
②用窗函数法设计FIR滤波器完成同步解调，观察滤波器频率响应曲线；（要求：
分别使用矩形窗和布莱克曼窗，进行比较）；
③输出解调信号的波形和频谱，观察现象，给出理论解释；
④播放解调音乐信号，比较不同滤波器下的声音，解释现象。
理论基础:
有一条信号f（t）恢复原始信号g（t）的过程称为解调。
这里，cos（w0*t）信号是接收端的本地载波信号，它与发送端的载波同频同相。
f（t）与cos（w0*t）相乘的结果使频谱F（W）向左、右分别移动+w0、-w0（并乘以系数1/2），得到
g0（t）=1/2*g（t）+1/2*g（t）*cos（w0*t）
G0（w）=1/2*G（w）+1/4*[G（w-2*w0）+G（w+2*w0）]
再利用一个低通滤波器，滤除在频率为2*w0附近的分量，即可取出g（t）,完成解调。
【Matlab程序如下】：
%理想低通滤波器冲击响应函数
functionhd=ideal（N,wc）
forn=0:
N-1
ifn==（N-1）/2
hd（n+1）=wc/pi;
elsehd（n+1）=sin（wc*（n-（N-1）/2））/（pi*（n-（N-1）/2））;
end
end
（将上述程序保存为ideal.m，但不能运行。
）
Clc
closeall
clearall
[y,fs]=wavread（'一生有你'）;
y1=y（:
1）;
N0=length（y）;
n=0:
（N0-1）;
w=2*n/（N0-1）;
t=n/fs;
f1=fft（y1）;
%调制频率滤波
y2=cos（n*0.4*pi）;%调制频率为0.4*pi;
y3=y1.*y2';%调制后的信号;
y4=y3.*y2';%调制解调后的信号;
%对调制频率巴特沃斯IIR滤波器的设计
Wp=0.2;
Ws=0.3;
Rp=1;
Rs=20;
W=0:
0.001*pi:
1*pi;
[N,Wc]=buttord（Wp,Ws,Rp,Rs）;
[b,a]=butter（N,Wc）;
[H,W]=freqz（b,a,W）;
plot（W/pi,abs（H））
title（'巴特沃斯滤波器的频域图'）
xlabel（'w'）
ylabel（'H（w）'）
y5=filter（b,a,y4）*2;%对调制解调后的信号进行滤波;
f3=fft（y3）;%对调制信号做fft变换;
f4=fft（y4）;%对调制解调后的信号做fft变换;
f5=fft（y5）;%对调制解调滤波后的信号做fft变换;
figure
subplot（4,1,1）,plot（w,abs（f1）*2/N0）%调制前的原信号频域图；
title（'原信号频域图'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（4,1,2）,plot（w,abs（f3）*2/N0）%高调制后的信号频域图；
title（'调制后的频域图'）
xlabel（'w'）
ylabel（'Y3（w）'）
subplot（4,1,3）,plot（w,abs（f4）*2/N0）
title（'调制解调后的频域图'）
xlabel（'w'）
ylabel（'Y4（w）'）
subplot（4,1,4）,plot（w,abs（f5）*2/N0）%解调滤波后的信号频域图；
title（'解调滤波后的频域图'）
xlabel（'w'）
ylabel（'Y5（w）'）
%用窗函数设计FIR滤波器
N=33;
Wc1=0.4*pi;
hd=ideal（N,Wc1）;
w1=boxcar（N）;%矩形窗
w2=blackman（N）;%布莱克曼窗
h1=hd.*w1';
h2=hd.*w2';
%求滤波器频率响应
M=1024;
W=2/M*[0:
M-1];
fh1=fft（h1,M）;
db1=-20*log10（abs（fh1
（1）./（abs（fh1）+eps）））;%幅度的分贝表示
fh2=fft（h2,M）;
db2=-20*log10（abs（fh2
（1）./（abs（fh2）+eps）））;%幅度的分贝表示
figure
subplot（2,2,1）;plot（W,abs（fh1））
title（'矩形窗幅度'）
subplot（2,2,2）;plot（W,db1）
title（'矩形窗分贝'）
subplot（2,2,3）;plot（W,abs（fh2））
title（'布莱克曼窗幅度'）
subplot（2,2,4）;plot（W,db2）
title（'布莱克曼窗分贝'）
y6=conv（y4,h1）;%作卷积
f6=fft（y6）;
y7=conv（y4,h2）;
f7=fft（y7）;
M1=length（y4）+length（h1）-1;
W1=2/M1*[0:
M1-1];
M2=length（y4）+length（h2）-1;
W2=2/M2*[0:
M2-1];
figure
subplot（3,1,1）,plot（w,abs（f1）*2/（length（y）））
title（'原信号频域图'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（3,1,2）,plot（W1,abs（f6）*2/M1）
title（'用矩形窗滤波后的频域图'）
xlabel（'w'）
ylabel（'Y6（w）'）
subplot（3,1,3）,plot（W2,abs（f7）*2/M2）
title（'用布莱克曼窗滤波后的频域图'）
xlabel（'w'）
ylabel（'Y7（w）'）
wavplay（y1,fs）
wavplay（y6,fs）
wavplay（y7,fs）
程序运行如下：
分析1：
通过对比解调前后的音乐信号时域图和频谱图可知，解调后的音乐信号时域和频谱幅度变为原来的1/2，若想滤波后可以完全恢复原始音乐信号，可以在传输函数添加系数2即可。
分析2：
对FIR窗函数滤波器的设计参数进行选定：
矩形窗（boxcar）：
滤波器的截止频率wc=0.5pi，滤波器的阶数N=33，阻带衰减为25db；布莱克曼窗（blackman）：
滤波器的截止频率wc=0.5pi，滤波器的阶数N=33，阻带衰减为80db。
比较采用矩形窗和布莱克曼窗的频率特性图可以看出：
最小阻带衰减只由窗形状决定，而不受阶数N的影响；而过渡带的宽度窗的形状有关；同时，布莱克曼窗的过滤带宽、旁瓣峰值和主瓣宽度均大于矩形窗的过滤带宽、旁瓣峰值和主瓣宽度。
分析3：
下图为经过加矩形窗和布莱克曼窗的滤波器滤波后的音乐信号的时域图和频谱图。
比较两种窗函数滤波下的频谱图可知加矩形窗的滤波器的滤波效果明显没有加布莱克曼窗的滤波器滤波效果好，原因如分析2所示：
布莱克曼窗的最小阻带衰减在同样阶数的前提条件下要比矩形窗的最小阻带衰减大，所以滤波更彻底，恢复的原信号也就会比矩形窗好。
若想提高矩形窗的滤波效果可采取增加阶数的方法，但这种情况下会增加不必要的开支，所以在选定滤波器的时候要权衡各方面的条件，力争达到一个各方都满意的结果。
分析4：
播放IIR滤波器和FIR滤波器滤波后的音乐信号，我们可以听到音乐信号的音调变低了，这是因为高频部分被滤掉了。
同时还可以听到IIR滤波器滤波后的音乐信号没有FIR滤波器滤波后的音乐信号效果好，这是因为FIR滤波器是线性相位的，而IIR滤波器是非线性相位；矩形窗没有布莱克曼窗的滤波效果好。
5、音乐信号的滤波去噪
①给原始音乐信号叠加幅度为0.05，频率为3kHz、5kHz、8kHz的三余弦混合噪声，观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；
②给原始音乐信号叠加幅度为0.5的随机白噪声（可用rand语句产生），观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；
③根据步骤①、②观察到的频谱，选择合适指标设计滤波器进行滤波去噪，观察去噪后信号音谱和频谱，并播放音乐，解释现象。
理论基础：
Rand函数介绍：
rand函数产生由在（0，1）之间均匀分布的随机数组组成的数组。
Y=rand（n）返回一个nxn的随机矩阵。
如果n不是数量，则返回错误信息。
【Matlab程序如下】：
Clc
Closeall
clearall
[y,fs]=wavread（'一生有你'）;
y1=y（:
1）;
f1=fft（y1）;
N0=length（y）;
n=0:
（N0-1）;
w=2*n/N0;
t=n/fs;
y2=0.05*cos（2*pi*3000*t）+0.05*cos（2*pi*5000*t）+0.05*cos（2*pi*8000*t）;%三余弦函数
f2=fft（y2）;
y3=y1+y2';
f3=fft（y3）;
figure
subplot（3,2,1）,plot（t（1:
1000）,y1（1:
1000））
title（'原信号时域图'）
xlabel（'t'）
ylabel（'y1'）
subplot（3,2,2）,plot（w,abs（f1）*2/N0）
title（'原信号频谱'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（3,2,3）,plot（t（1:
200）,y2（1:
200））
title（'三余弦信号时域图'）
xlabel（'t'）
ylabel（'y2'）
subplot（3,2,4）,plot（w,abs（f2）*2/N0）
title（'三余弦信号频谱'）
xlabel（'w'）
ylabel（'Y2（w）'）
subplot（3,2,5）,plot（t（1:
1000）,y3（1:
1000））
title（'加噪后信号的时域'）
xlabel（'t'）
ylabel（'y2'）
subplot（3,2,6）,plot（w,abs（f3）*2/N0）
title（'加噪后信号的的频谱'）
xlabel（'w'）
ylabel（'Y3（w）'）
y4=rand（N0,1）-0.5;%白噪声的产生算法
y5=y1+y4;
f4=fft（y4）;
f5=fft（y5）;
figure
subplot（2,2,1）,plot（t（1:
1000）,y4（1:
1000））
title（'随机白噪声时域图'）
xlabel（'t'）
ylabel（'y4'）
subplot（2,2,3）,plot（t（1:
1000）,y5（1:
1000））
title（'带有随机白噪声的信号时域'）
xlabel（'t'）
ylabel（'y5'）
subplot（2,2,2）,plot（w,abs（f4）*2/N0）
title（'随机白噪声频谱'）
xlabel（'w'）
ylabel（'Y4（w）'）
subplot（2,2,4）,plot（w,abs（f5）*2/N0）
title（'带有随机白噪声的信号频谱'）
xlabel（'w'）
ylabel（'Y5（w）'）
%巴特沃斯滤波器的设计
Wp=0.08;
Ws=0.12;%不宜过狭窄不然会出错
Rp=1;
Rs=30;%指标输入
[N,Wc]=buttord（Wp,Ws,Rp,Rs）;
[b,a]=butter（N,Wc）;
[H,W]=freqz（b,a）;
figure
plot（W/pi,abs（H））
title（'低巴特沃斯滤波器的频域图'）
xlabel（'w'）
ylabel（'H（w）'）
y6=filter（b,a,y3）;
y7=filter（b,a,y5）;
f6=fft（y6）;
f7=fft（y7）;
figure
subplot（3,1,1）,plot（w,abs（f1）*2/N0）
title（'原信号频谱'）
xlabel（'w'）
ylabel（'Y1（w）'）
subplot（3,1,2）,plot（
升级会员 