数字信号处理课程设计基于 matlab 的音乐信号处理和分析.docx
《数字信号处理课程设计基于 matlab 的音乐信号处理和分析.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计基于 matlab 的音乐信号处理和分析.docx（36页珍藏版）》请在冰豆网上搜索。
数字信号处理课程设计基于matlab的音乐信号处理和分析
《数字信号处理》课程设计
设计题目：
基于MATLAB的音乐信号处理和分析
一、课程设计的目的
本课程设计通过对音乐信号的采样、抽取、调制、解调等多种处理过程的理论分析和MATLAB实现，使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；使学生掌握的基本理论和分析方法知识得到进一步扩展；使学生能有效地将理论和实际紧密结合；增强学生软件编程实现能力和解决实际问题的能力。
二、课程设计基本要求
1学会MATLAB的使用，掌握MATLAB的基本编程语句。
2掌握在Windows环境下音乐信号采集的方法。
3掌握数字信号处理的基本概念、基本理论和基本方法。
4掌握MATLAB设计FIR和IIR数字滤波器的方法。
5掌握使用MATLAB处理数字信号、进行频谱分析、设计数字滤波器的编程方法。
三、课程设计内容
1、音乐信号的音谱和频谱观察
使用windows下的录音机录制一段音乐信号或采用其它软件截取一段音乐信号（要求：
时间不超过5s、文件格式为wav文件）
①使用wavread语句读取音乐信号，获取抽样率；（注意：
读取的信号是双声道信号，即为双列向量，需要分列处理）；
②输出音乐信号的波形和频谱，观察现象；
③使用sound语句播放音乐信号，注意不同抽样率下的音调变化，解释现象。
Wavread格式说明：
[w,fs,b]=wavread（‘语音信号’），采样值放在向量w中，fs表示采样频率（hz），b表示采样位数。
上机程序：
[y,fs,bit]=wavread（'Ido片段'）%读取音乐片段，fs是采样率
size（y）%求矩阵的行数和列数
y1=y（:
1）;%对信号进行分列处理
n1=length（y1）;%取y的长度
t1=（0:
n1-1）/fs;%设置波形图横坐标
figure
subplot（2,1,1）;
plot（t1,y1）;%画出时域波形图
ylabel（'幅值'）;
xlabel（'时间（s）'）;
title（'信号波形'）;
subplot（2,1,2）;
Y1=fft（y1）;
w1=2/n1*（0:
n1-1）;%设置角频率
plot（w1,abs（Y1））;%画频谱图
title（'信号频谱'）;
xlabel（'数字角频率'）;
ylabel（'幅度'）;
gridon;
sound（y,fs）;
实验结果：
1、通过观察频谱知，选取音乐信号的频谱集中在0~0.7*pi之间，抽样点数fs=44100;
2、当采样频率问原来0.5（0.5*fs）倍时：
音乐片段音调变得非常低沉，无法辨认原声，播放时间变长；抽样频率减小，抽样点数不变时，其分辨力增大，记录长度变长，声音失真。
3、当采样频率问原来2（2*fs）倍时：
音乐片段音调变得尖而细，语速变快，播放时间变短；抽样频率增加，抽样点数不变时，其分辨力下降，记录长度变短，声音失真。
2、音乐信号的抽取（减抽样）
①观察音乐信号频率上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；
②输出减抽样音乐信号的波形和频谱，观察现象，给出理论解释；
③播放减抽样音乐信号，注意抽样率的改变，比较不同抽取间隔下的声音，解释现象。
上机程序：
[y,fs,bit]=wavread（'Ido片段'）
y1=y（:
1）;
n1=length（y1）;
tn1=（0:
n1-1）/fs;
figure
subplot（2,1,1）;
plot（tn1,y1）;
ylabel（'幅度'）;
xlabel（'时间（s）'）;
title（'原信号波形'）;
wn1=2/n1*[0:
n1-1];
Y1=fft（y1）;
subplot（2,1,2）;
plot（wn1,abs（Y1））;
title（'原信号频谱'）;
xlabel（'数字角频率w'）;
ylabel（'幅度'）;
gridon;
D=2;%设置抽样间隔
y2=y1（1:
D:
n1）;%减抽样
n2=length（y2）;%减抽样后信号长度
t2=（0:
n2-1）/fs;%设置横坐标
figure
subplot（2,1,1）;
plot（t2,y2）;%绘制减抽样信号波形图
ylabel（'幅度'）;
xlabel（'时间（s）'）;
title（'2:
1减抽样信号波形'）;
Y2=fft（y2）;%对y2进行n2点fft谱分析
w2=2/n2*[0:
n2-1];
subplot（2,1,2）;
plot（w2,abs（Y2））;%绘制减抽样信号频谱图
title（'2:
1减抽样信号频谱'）;
xlabel（'数字角频率w'）;
ylabel（'幅度'）;
gridon;
sound（y2,fs/D）;
实验结果与分析：
1、程序中指标D表示抽样间隔，其值越大，相邻两抽样点之间的距离越远，抽样后漏掉的信息越多，相应的时域信号长度越短；
2、抽样间隔D=1.1时的信号波形及频谱图，抽样频率
大于信号最高频率
的两倍，
满足抽样定理，不会发生混叠。
抽样间隔D越大，抽样率fs越小，抽样后时域信号长度越短
3、抽样间隔D=2时的信号波形及频谱图
抽样间隔D=2的信号波形及频谱图，抽样频率
小于信号最高频率
的两倍，即
<2
，不满足抽样定理，其频谱图发生混叠，且D越大，混叠越严重，高频成分增加越多，音乐片段音调听起来很沙哑，音调变得很高。
4、抽样间隔D=5的信号波形及频谱图，抽样频率
小于信号最高频率
的两倍，即
<2
，不满足抽样定理，其频谱图发生混叠，音乐片段音调听起来很尖锐。
3、音乐信号的AM调制
①观察音乐信号频率上限，选择适当调制频率对信号进行调制（给出高、低两种调制频率）；
②输出调制信号的波形和频谱，观察现象，给出理论解释；
③播放调制音乐信号，注意不同调制频率下的声音，解释现象。
上机程序：
[y,fs,bit]=wavread（'Ido片段'）
y1=y（:
1）;
n1=length（y1）;
n3=0:
（n1-1）;
b1=cos（0.75*pi*n3）;%设置调制信号
c1=b1'.*y1;%对原信号进行调制
lc1=length（c1）;
t=（0:
lc1-1）/fs;
figure%用载波对信号进行调制，并对其做fft变换
subplot（2,1,1）%获取频谱，从图中可以观察到，调制后的
plot（t,c1）;%信号频谱发生搬移
xlabel（'时间（s）'）;
ylabel（'幅度'）;
title（'调制后信号'）;
w1=2/lc1*[0:
lc1-1];%设置角频率W
C1=fft（c1）;
subplot（2,1,2）
plot（w1,abs（C1））;
xlabel（'数字角频率w'）;
ylabel（'幅度'）;
title（'调制后信号的频谱（高频率调制）'）;
gridon;
%sound（c1,fs）;
实验结果与分析：
信号的调制过程就是将信号频谱搬移到任何所需的较高频率范围。
调制的实质是把各种信号的频谱搬移，使它们互不重叠的占据不同的频率范围，也即信号分别托付于不同频率的载波上。
具体的调制原理推导在此不再叙述，仅将结论列出：
f（t）=g（t）*cos（w0*t）。
由此将信号g（t）的频谱搬移到（2*n+1）w附近，同时音乐信号的频谱幅度变为原来的1/2。
如果信号的最高频谱wh超过了ws/2，则各周期延括分量产生频谱的交叠，称为是频谱的混叠现象。
根据奈奎斯特定律可知，若希望频谱不会发生混叠，则fs>=fh。
调制后信号波形及频谱（低频率调制）b1=cos（0.25*pi*n3）;
调制后信号波形及频谱（高频率调制）b1=cos（0.75*pi*n3）;
1、由频谱图知信号频率上限约在0.7*pi处，取高调制频率为0.75*pi，取高调制频率为0.25*pi；观察调制后的频谱知，高频率调制后频谱混叠，而低频率调制后频谱未发生混叠；
2、用载波对信号进行调制，其实质过程为在时域用载波函数b1=cos（0.25*pi*n3）乘上原信号函数，在频域既表现为信号频谱的搬移，且载波函数角频率越大，调制后信号频谱搬移越大；
3、调制后信号高频成分增加，且调制频率越高，高频成分越大，声音变得越低沉，失真程度越大。
4、AM调制音乐信号的同步解调
①设计巴特沃斯IIR滤波器完成同步解调；观察滤波器频率响应曲线；
②用窗函数法设计FIR滤波器完成同步解调，观察滤波器频率响应曲线；（要求：
分别使用矩形窗和布莱克曼窗，进行比较）；
③输出解调信号的波形和频谱，观察现象，给出理论解释；
④播放解调音乐信号，比较不同滤波器下的声音，解释现象。
上机程序：
clearall;cla;closeall
[a,fs,bit]=wavread（'Ido片段'）;
y1=a（:
1）;%去单列数据进行分析
f1=fft（y1）;
n=length（f1）;
tn=（0:
n-1）/fs;
w=2/n*[0:
n-1];
%sound（y1,fs）;
figure
（1）
subplot（2,2,1）;plot（tn,y1）;%绘制音乐信号时域波形图
gridon;
title（'原信号音频'）;
xlabel（'时间'）;
ylabel（'幅度'）;
subplot（2,2,2）;plot（w,abs（f1））;%绘制音乐信号频谱图
gridon;
title（'原信号频谱'）;
xlabel（'频率/pi'）;
ylabel（'幅度'）;
t=[0:
n-1];
y2=cos（pi*1/2*t）;%载波函数
y3=y1.*y2';%信号调制
ty3=（0:
length（y3）-1）/fs;
subplot（2,2,3）;plot（ty3,y3）;%绘制调制后信号波形图
gridon;
title（'AM调制音频信号'）;
xlabel（'时间'）;
ylabel（'幅度'）;
f3=fft（y3）;
n2=length（f3）;
w2=2/n2*[0:
n2-1];
subplot（2,2,4）;plot（w2,abs（f3））;%绘制调制后信号频谱图
gridon;
title（'AM调制频谱'）;
xlabel（'频率/pi'）;
ylabel（'幅度'）;
%解调后信号
y4=y3.*y2'%信号解调
ly4=length（y4）;
ty4=（0:
（ly4-1））/fs;
wy4=2/ly4*[0:
ly4-1];
Y4=fft（y4）;
figure
（2）
subplot（3,1,1）;plot（ty4,y4）;
gridon;
title（'解调后的波形'）;
xlabel（'t'）;
ylabel（'幅度'）;
subplot（3,1,2）;plot（wy4,abs（Y4））;
title（'解调后波形频谱'）;
xlabel（'w/pi'）;
ylabel（'幅度'）;
gridon;
%设计巴特沃斯滤波器进行滤波去噪
[N1,wc1]=buttord（0.05,0.50,1,15）;%确定低通滤波器的阶数和截止频率；
[b,a]=butter（N1,wc1）;%确定低通滤波器分子分母系数
[H,W]=freqz（b,a）;
subplot（3,1,3）;
plot（W,abs（H））;
title（'IIR滤波器频谱'）;
xlabel（'w/pi'）;
ylabel（'幅度'）;
gridon;
m=filter（b,a,y4）;
sound（m,fs）;
lm=length（m）;%滤波后信号长度
tm=（0:
lm-1）/fs;%设置横坐标
wm=2/lm*[0:
lm-1];
M=fft（m）;
figure（3）
subplot（2,1,1）;plot（tm,m）;
gridon;
title（'IIR滤波器滤波后波形'）;
xlabel（'t'）;
ylabel（'幅度'）;
subplot（2,1,2）;plot（wm,abs（M））;
title（'IIR滤波器滤波后频谱'）;
xlabel（'w/pi'）;
ylabel（'幅度'）;
%矩形窗和布莱克曼窗
N=33;wc=0.3*pi;%基于经验的指标，其中N为理想低通滤波器阶数，wc为截止频率
hd=ideal（N,wc）;%调用理想低通滤波器函数
w1=boxcar（N）;%产生各种窗函数
w2=blackman（N）;
h1=hd.*w1';%加窗设计各种FIR滤波器
h2=hd.*w2';
th1=（0:
32）/fs;
th2=（0:
32）/fs;
M=21184;
fh1=fft（h1,M）;%矩形窗频谱函数
w=2/M*[0:
M-1];
fh2=fft（h2,M）;%布莱克曼窗频谱函数
figure（4）
subplot（2,2,1）;plot（th1,h1）
title（'矩形窗时域'）;
subplot（2,2,2）;plot（w,abs（fh1））;
title（'矩形窗频域'）;
subplot（2,2,3）;plot（th2,h2）;
title（'布莱克曼窗时域'）;
subplot（2,2,4）;plot（w,abs（fh2））;
title（'布莱克曼窗频域'）
%滤波处理
y6=conv（h1,y4）;%用矩形窗对调制后信号进行滤波
f6=fft（y6）;
n4=length（f6）;
ty6=（0:
n4-1）/fs;
w3=2/n4*[0:
n4-1];
%sound（y6,fs）;
figure（5）
subplot（2,2,1）;plot（ty6,y6）;
title（'矩形窗滤波后音频'）
subplot（2,2,2）;plot（w3,abs（f6））;
title（'矩形窗滤波后频谱'）
y7=conv（h2,y4）;%用布莱克曼窗对调制后的信号进行滤波
f7=fft（y7）;
n5=length（f7）;
ty7=（0:
n5-1）/fs;
w4=2/n5*[0:
n5-1];
%sound（y7,fs）;
subplot（2,2,3）;plot（ty7,y7）;
title（'布莱克曼窗滤波后音频'）
subplot（2,2,4）;plot（w4,abs（f7））;
title（'布莱克曼窗滤波后频谱'）
实验结果与分析：
由调制信号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）,完成解调。
窗函数法的设计步骤：
a、给定所要求的频率响应函数Hd（w）；
b、求出hd（n）；
c、有过渡带宽及阻带最小衰减的要求，选定窗w（n）的形状以及N的大小，一般N要通过几次试探才能确定；
d、求得所设计的FIR滤波器的单位抽样响应：
h（n）=hd（n）*w（n）；
e、求H（w）,检验是否符合设计要求，如不满足，则需重新设计。
函数格式说明：
[n,wn]=buttord（wp,ws,rp,rs）：
n为滤波器的阶数，wn为其截止频率，wp为通带截止频率，ws为阻带截止频率，rp为通带最大衰减,，rs为阻带最小衰减。
[a,b]=butter（n,rp）：
a，b为所设计的滤波器传输函数分子和分母的系数。
[h,w]=freqz（a,b）:
求滤波器的频率响应h，w为频域坐标。
1、由调制后的信号频谱知调制过程实质为信号频域的搬移，而解调则为调制的逆过程；
根据解调后信号的频谱，设置巴特沃斯IIR低通滤波器的通带截止频率为0.5，阻带截止率为0.50，阻带和通带衰减分别为1db和15db;
2、由滤波器滤波后波形频谱知，与原信号频谱相比原信号频谱中的一些高频成分也被滤除，但总体来说滤波器所选指标还是能够较好地还原信号信息，通过播放滤波后音频也可以说明这一点；
3、对于窗函数的设计，其指标则是基于频谱图和过往经验的综合，为了便于比较两种窗函数滤波的效果，两种窗函数选用相同的指标；实验中设置窗函数的阶数N=33，其截止频率Wc=0.3*pi（比巴特沃斯IIR低通滤波器截止频率略大）；通过比较两窗函数的频谱和滤波后的频谱知，矩形窗通带与阻带的过渡部分比较陡，对于高频信息滤除的比较彻底；而布莱克曼窗的过渡部分则相对比较平缓，相应的其对于高频信息的滤除就没有矩形窗那么彻底（相对而言）；比较采用矩形窗和布莱克曼窗的频率特性图可以看出：
最小阻带衰减只由窗形状决定，而不受阶数N的影响；而过渡带的宽度窗的形状有关；同时，布莱克曼窗的过滤带宽、旁瓣峰值和主瓣宽度均大于矩形窗的过滤带宽、旁瓣峰值和主瓣宽度；分别播放两种窗函数滤波后的声音，通过布莱克曼窗的声音比通过矩形窗的声音稍微响亮一些（不是很明显），这也印证了上面对于两种窗函数的分析。
5、音乐信号的滤波去噪
①给原始音乐信号叠加幅度为0.05，频率为3kHz、5kHz、8kHz的三余弦混合噪声，观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；
②给原始音乐信号叠加幅度为0.5的随机白噪声（可用rand语句产生），观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；
③根据步骤①、②观察到的频谱，选择合适指标设计滤波器进行滤波去噪，观察去噪后信号音谱和频谱，并播放音乐，解释现象。
源程序:
clc;clear;close;
[y,fs,bit]=wavread（'Ido片段'）;
y0=y（:
1）;
l=length（y0）;
%加三余弦混合噪声
t0=（0:
l-1）/fs;
d0=[0.05*cos（2*pi*3000*t0）]';
t1=（0:
l-1）/fs;
d1=[0.05*cos（2*pi*5000*t1）]';
t2=（0:
l-1）/fs;
d2=[0.05*cos（2*pi*8000*t2）]';
noise=d2+d1+d0;
y1=y0+noise;
%sound（y1,fs）;
a=length（noise）;%绘制三余弦噪声音频图
wa=2/a*[0:
a-1];
Noise=fft（noise）;
figure
（1）
subplot（2,2,3）;
plot（noise（1:
150））;
xlabel（'时间（s）'）
ylabel（'幅值'）
title（'三余弦信号音谱'）
subplot（2,2,1）;%绘制三余弦噪声频谱图
plot（wa,abs（Noise））;
gridon;
xlabel（'W'）
title（'噪声频谱'）
w0=2/l*[0:
l-1];%绘制加噪信号音频
Y1=fft（y1）;
subplot（2,2,4）
plot（w0,abs（Y1））;
gridon;
xlabel（'W'）
title（'加噪信号频谱'）
ly1=length（y1）;
ty1=（0:
ly1-1）/fs;
subplot（2,2,2）;plot（ty1,y1）;
xlabel（'时间（s）'）
ylabel（'幅值'）
title（'加噪信号音谱'）
m=rand（l,1）-0.5;%产生幅度为0.5的随机信号
lm=length（m）;
y2=m+y0;%将噪声信号与原声音信号叠加
%sound（y2,fs）;
wm=2/lm*[0:
lm-1];
M=fft（m）;
figure
（2）
subplot（2,2,3）;
plot（m（1:
150））
xlabel（'时间（s）'）
ylabel（'幅值'）
title（'白噪信号音谱'）
subplot（2,2,1）;
plot（wm,abs（M））;
gridon;
xlabel（'W'）
title（'噪声频谱'）
l=length（y2）;
ty2=（0:
l-1）/fs;
w=2/l*[0:
l-1];
Y2=fft（y2）;
subplot（2,2,4）
plot（w,abs（Y2））;
gridon;
xlabel（'W'）
title（'加噪信号频谱'）
subplot（2,2,2）;plot（ty2,y2）;
xlabel（'时间（s）'）
ylabel（'加噪信号幅值'）
title（'加噪信号音谱'）;
%设计滤波器进行滤波去噪
[N1,wc1]=buttord（0.04,0.13,1,24）;%确定低通滤波器的阶数和截止频率；
%[N1,wc1]=buttord（0.04,0.20,1,18）;（针对随机白噪声的指标）
[b,a]=butter（N1,wc1）;%确定低通滤波器分子分母系数
[H,W]=freqz（b,a）;
figure（3）;
subplot（3,1,1）;
plot（W,abs（H））;%滤波器频谱
xlabel（'w/pi'）
ylabel（'幅度k'）
title（'IIR滤波器滤波波形'）;
gridon;
m=filter（b,a,y1）;%用滤波器滤除三余弦噪声
%sound（m,fs）;
lm=length（m）;%滤波后信号长度
tm=（0:
lm-1）/fs;%设置横坐标
subplot（3,1,2）;
plot（tm,m）;%绘制滤波后的波形
xlabel（'t（s）'）
ylabel（'信号幅值'）
title（'去噪后信号波形'）;
k=fft（m）;%滤波后的波形做离散傅里叶变换
w=2*[0:
length（k）-1]/length（k）;
subplot（3,1,3）
plot（w,abs（k））;
xlabel（'w/pi'）
ylabel（'幅度k'）
title（'IIR滤波器滤波后信号频谱'）;
实验结果与分析：
信号的去噪即根据信号噪声频谱的分布范围，从而采用不同的滤波器对其进行
滤波。
常用的滤波器有IIR巴特沃斯低通滤波器，FIR窗函数滤波器。
Rand函数介绍：
rand函数产生由在（0，1）之间均匀分布的随机数组组成的数组。
Y=rand（n）返回一个nxn的随机矩阵。
如果n不是数量，则返回错误信息。
Y=rand（m,n）或Y=rand（[mn]）返回一个mxn的随机矩阵。
Y=rand（m,n,p,...）或Y=rand（[mnp...]）产生随机数组。
Y=rand（size（A））返回一个和A有相同尺寸的随机矩阵。
1、rand（3）*-2 ：
rand（3）是一个3*3的随机矩阵（数值范围在0~1之间）然后就是每个数乘上-2；2、
升级会员 