DSP课程设计.docx
《DSP课程设计.docx》由会员分享，可在线阅读，更多相关《DSP课程设计.docx（24页珍藏版）》请在冰豆网上搜索。
DSP课程设计
《数字信号处理》
课程设计报告
题目:
基于MATLAB的音乐信号处理和分析
专业：
测控技术与仪器
年级：
2011级
姓名：
学号：
基于MATLAB的音乐信号处理和分析
一、读取音乐信号并观察波形与频谱
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））';
sound（w,fs）
figure;
subplot（2,1,1）;plot（wav）;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;f=wwav/2*fs;
subplot（2,1,2）;
plot（wwav,abs（fwav（1:
lwav）））;
1，以二倍的抽样率听声音信号时，音乐播放的特别快，像被压缩了，播放的时间比原信号短。
2，以二分之一的抽样率听声音信号时，音乐播放的特别慢，像被拉长了，播放的时间比原信号长。
3，原信号频谱截止频率为0.5*pi。
二、
（1）音乐信号的5倍减抽样
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））';
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
j=0;d1=5;
fori=1:
d1:
lwav
j=j+1;
dwav1（j）=wav（i）;
end
figure;
subplot（2,1,1）;
plot（dwav1）;
fdwav1=fft（dwav1）;
sound（dwav1,fs/d1）
lwav1=round（length（fdwav1）/2）;
nwav=[0:
lwav1-1];
wwav=nwav/（lwav1）;f=wwav/2*fs;
subplot（2,1,2）;
plot（wwav,abs（fdwav1（1:
lwav1）））;%减抽样音乐频域
（2）音乐信号的10倍减抽样（有混频）
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））';
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
j=0;d1=10;
fori=1:
d1:
lwav
j=j+1;
dwav1（j）=wav（i）;
end
figure;
subplot（2,1,1）;
plot（dwav1）;
fdwav1=fft（dwav1）;
sound（dwav1,fs/d1）
lwav1=round（length（fdwav1）/2）;
nwav=[0:
lwav1-1];
wwav=nwav/（lwav1）;f=wwav/2*fs;
subplot（2,1,2）;
plot（wwav,abs（fdwav1（1:
lwav1）））;
1，原信号频谱截止频率为0.5*pi。
2，减抽样后的音乐信号听起来变得尖锐，有失真。
3，抽样率随着抽样间隔的增大而逐渐变小，声音越来越失真，音调变得急促，而尖锐，信号产生混叠
3、音乐信号的AM调制
（1）调制频率为pi/2
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））';
lam=length（wav）;
amw=pi/2;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav;
sound（amwav,fs）
figure;
subplot（2,1,1）;plot（amwav）;
title（'爸爸去哪儿调制后'）
famwav=fft（amwav）;
lwav=round（length（famwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
subplot（2,1,2）;
plot（wwav,abs（famwav（1:
lwav）））;
（2）调制频率为2*pi/3
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））';
lam=length（wav）;
amw=2*pi/3;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav;
sound（amwav,fs）
figure;
subplot（2,1,1）;plot（amwav）;
title（'爸爸去哪儿调制后'）
famwav=fft（amwav）;
lwav=round（length（famwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
subplot（2,1,2）;
plot（wwav,abs（famwav（1:
lwav）））;
1，原信号的调制相当于频谱搬移，左移一个右移一个，当调制频率（余弦频率）小于0.4pi或大于0.6pi时就会产生混叠或丢失一部分信息。
2，当余弦点数取得少时，余弦频谱会产生泄漏。
3，当调制频率较高时（发生混叠），声音响度低，几乎只能听见兹兹的声音，信号几乎完全失真，当调制频率较低时（未发生混叠），声音很尖锐，响度较大，稍微能听出一点调子，但也有兹兹的声音。
四、AM调制音乐信号的同步解调
1.信号的解调
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
lam=length（wav）;
amw=pi/2;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav';
pmwav=amwav.*cwav';%音乐信号解调
figure;
subplot（2,1,1）;plot（pmwav）;
title（'解调时域'）
fpmwav=fft（pmwav）;
subplot（2,1,2）;
plot（wwav,abs（fpmwav（1:
lwav）））;
title（'解调频域'）
sound（pmwav,fs）;
2.巴特沃斯滤波
[w,fs,bit]=wavread（'爸爸去哪儿.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
lam=length（wav）;
amw=pi/2;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav';
pmwav=amwav.*cwav';%音乐信号解调
fpmwav=fft（pmwav）;
T=1;
[N,wc]=buttord（0.05,0.5,1,15,'s'）;
[B,A]=butter（N,wc,'s'）;%求模拟滤波器系统函数的系数
[C,D]=bilinear（B,A,1/T）;%求数字滤波器系统函数的系数
W=0:
0.001*pi:
0.5*pi;
H=freqs（B,A,W）;%模拟滤波器的频率响应
Hd=freqz（C,D,W）;%数字滤波器频率响应
figure;
Hd_db=-20*log（abs（Hd
（1）./（abs（Hd）+eps）））;
subplot（3,1,1）;plot（W/pi,abs（Hd））;
title（'数字巴特沃斯滤波频响图'）
gridon;
y=filter（C,D,pmwav）;%对信号滤波
fy=fft（y）;
subplot（3,1,2）;plot（abs（y））;
title（'滤波后音频'）
subplot（3,1,3）;plot（wwav,abs（fy（1:
lwav）））;
title（'滤波后频谱'）
sound（y,fs）;
3.矩形窗滤波
[w,fs,bit]=wavread（'爸爸.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
lam=length（wav）;
amw=pi/2;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav';
pmwav=amwav.*cwav';%音乐信号解调
fpmwav=fft（pmwav）;
N=39,wc=0.35*pi;%参数设置
M=512;
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;%调用理想低通滤波器函数
w2=boxcar（N）;%设置矩形窗
h2=hd.*w2';%求数字滤波器频率响应
lx=round（length（h2）/2）;
wx1=（2/M）*[0:
M-1];
fh2=fft（h2,M）;
db1=-20*log10（abs（fh2
（1）./（abs（fh2）+eps）））;
figure;
subplot（4,1,1）;plot（h2）;
title（'矩形窗滤波特性'）
gridon;
subplot（4,1,2）;plot（wx1,abs（fh2））;
title（'矩形窗频谱'）
gridon;
y0=pmwav;
yy2=conv（y0,h2）;
fyy2=fft（yy2）;
l2=round（length（fyy2）/2）;
wx2=（0:
l2-1）/l2;
subplot（4,1,3）;plot（wx2,abs（yy2（1:
l2）））;
title（'矩形窗滤波后音谱'）
gridon;
subplot（4,1,4）;plot（wx2,abs（fyy2（1:
l2）））;
title（'矩形窗滤波后频谱'）
gridon;
sound（yy2,fs）;
4、布莱克曼窗滤波
[w,fs,bit]=wavread（'爸爸.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
lam=length（wav）;
amw=pi/2;
cwav=cos（amw*[0:
lam-1]）;
amwav=wav.*cwav';
pmwav=amwav.*cwav';%音乐信号解调
fpmwav=fft（pmwav）;
N=43,wc=0.3*pi;%参数设置
M=512;
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;%调用理想低通滤波器函数
w2=blackman（N）;%设置矩形窗
h2=hd.*w2';%求数字滤波器频率响应
fh2=fft（h2）;
lx=length（h2）;
wx1=（0:
lx-1）*2/lx;
db1=-20*log10（abs（fh2
（1）./（abs（fh2）+eps）））;
figure;
subplot（4,1,1）;plot（h2）;
title（'布莱克曼窗滤波特性'）
gridon;
subplot（4,1,2）;plot（wx1,abs（fh2））;
title（'布莱克曼窗频谱'）
gridon;
y0=pmwav;
yy2=conv（y0,h2）;
subplot（4,1,3）;plot（yy2）;
title（'布莱克曼窗滤波后音谱'）
gridon;
fyy2=fft（yy2）;
l2=round（length（fyy2）/2）;
wx2=（0:
l2-1）/l2;
subplot（4,1,4）;plot（wx2,abs（fyy2（1:
l2）））;
title（'布莱克曼窗滤波后频谱'）
gridon;
sound（yy2,fs）;
1，解调后信号频谱在高频和低频处均有一部分，且成对称分布，需要滤掉高频才可大致还原原信号。
2，原信号的截止频率为0.5pi，使用数字巴特沃斯滤波器滤波器滤波参数通带截止频率0.5pi，阻带开始频率0.5pi，阻带衰减15db。
滤波效果很好，基本还原了原信号。
3，使用窗函数滤波要根据过渡带宽算阶数N，选截止频率为0.5pi。
4，使用矩形窗滤波，矩形窗过渡带窄，但是阻带有波纹，高频部分有小部分未滤掉。
5，使用布莱克曼窗滤波，布莱克曼窗过渡带宽，但是阻带较好。
6，使用矩形窗和布莱克曼窗滤波，效果都行，基本都能还原原信号
5、音乐信号的滤波去噪
1.原始信号加入三余弦混合噪声
[w,fs,bit]=wavread（'爸爸.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
sf1=3000;sf2=5000;sf3=8000;
l=length（wav）;
T=2/fs;
t=0:
T:
（l-1）*T;
s1=0.05*（cos（2*pi*sf2*t）+cos（2*pi*sf2*t）+cos（2*pi*sf3*t））;%余弦噪声信号
wav_s1=wav'+s1;%加噪信号
sound（wav_s1,fs）
fwav_s1=fft（wav_s1）;
f_s1=fft（s1）;
figure;
subplot（2,2,1）;plot（s1）;
title（'余弦噪声'）
subplot（2,2,2）;plot（wav_s1）;
title（'加噪信号'）
subplot（2,2,3）;
plot（wwav,abs（f_s1（1:
lwav）））;
subplot（2,2,4）;
plot（wwav,abs（fwav_s1（1:
lwav）））;
2.用窗函数对混合噪声滤波
[w,fs,bit]=wavread（'爸爸.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
sf1=3000;sf2=5000;sf3=8000;
l=length（wav）;
T=2/fs;
t=0:
T:
（l-1）*T;
s1=0.05*（cos（2*pi*sf2*t）+cos（2*pi*sf2*t）+cos（2*pi*sf3*t））;%余弦噪声信号
wav_s1=wav'+s1;%加噪信号
fwav_s1=fft（wav_s1）;
f_s1=fft（s1）;
N=45;wc=0.15*pi;%参数设置
M=512;
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;%调用理想低通滤波器函数
w1=hamming（N）;%海明窗
h1=hd.*w1';
fh1=fft（h1）;
lx=length（h1）;
wx1=[0:
lx-1]*2/lx;
db1=-20*log10（abs（fh1
（1）./（abs（fh1）+eps）））;
figure;
subplot（4,1,1）;plot（h1）;
title（'海明窗波形'）
gridon;
subplot（4,1,2）;plot（wx1,fh1）;
title（'海明窗频谱'）
gridon;
y0=wav_s1;
yy1=conv（y0,h1）;
subplot（4,1,3）;plot（yy1）;
title（'海明窗滤波后音频'）
gridon;
fyy1=fft（yy1）;
l2=round（length（fyy1）/2）;
wx2=（0:
l2-1）/l2;
subplot（4,1,4）;
plot（wx2,abs（fyy1（1:
l2）））;
title（'海明窗滤波后频谱'）
gridon;
sound（yy1,fs）
3.原始信号叠加随机白噪声并对其滤波
[w,fs,bit]=wavread（'爸爸.wav'）;
wav=（w（:
1））;
fwav=fft（wav）;
lwav=round（length（fwav）/2）;
nwav=[0:
lwav-1];
wwav=nwav/（lwav）;
l=length（wav）;
s2=（rand（l,1）-0.5）/5;%用rand语句设置幅度为0.5的随机噪声
wav_s2=wav+s2;
fwav_s2=fft（wav_s2）;
f_s2=fft（s2）;
figure;
subplot（2,2,1）;plot（s2）;
title（'随机白噪声'）
subplot（2,2,2）;plot（wav_s2）;
title（'加噪噪声'）
subplot（2,2,3）;plot（wwav,abs（f_s2（1:
lwav）））;
subplot（2,2,4）;plot（wwav,abs（fwav_s2（1:
lwav）））;
N=133;wc=0.15*pi;%参数设置
M=512;
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;%调用理想低通滤波器函数
w2=blackman（N）;%布莱克曼窗
h2=hd.*w2';
fh2=fft（h2）;
lx=length（h2）;
wx1=[0:
lx-1]*2/lx;
db1=-20*log10（abs（fh2
（1）./（abs（fh2）+eps）））;
figure;
subplot（4,1,1）;plot（wx1,abs（h2））;
title（'blackman窗波形'）
subplot（4,1,2）;plot（wx1,abs（fh2））;
title（'blackman窗频谱'）
gridon;
y0=wav_s2;
yy2=conv（y0,h2）;
subplot（4,1,3）;plot（yy2）;
title（'blackman窗滤波后音频'）
gridon;
fyy2=fft（yy2）;
l2=round（length（fyy2）/2）;
wx2=（0:
l2-1）/l2;
subplot（4,1,4）;plot（wx2,abs（fyy2（1:
l2）））;
title（'blackman窗滤波后频谱'）
gridon;
sound（yy2,fs）
1.三余弦信号的频谱为不同频率处得三根线，加噪声后的音乐信号频谱是在原信号频谱上加了三条不同频率的线。
加噪声后音乐信号能听到原有的音调，但里面有非常大的杂音，兹兹的噪声。
2.对加余弦噪声的信号进行滤波，用海明窗滤波，滤波后信号基本没有杂音，有稍微的失真。
3.对原信号加随机白噪声，白噪声均匀分布，对其用布莱克曼窗滤波，滤波后
原声还很清楚，但伴随有较大声的随机白噪声。
升级会员 