数字信号处理课程设计基于 matlab 的音乐信号处理和分析Word文档格式.docx

数字信号处理课程设计基于 matlab 的音乐信号处理和分析Word文档格式.docx

《数字信号处理课程设计基于 matlab 的音乐信号处理和分析Word文档格式.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计基于 matlab 的音乐信号处理和分析Word文档格式.docx（36页珍藏版）》请在冰豆网上搜索。

%画出时域波形图

ylabel（'

幅值'

）;

xlabel（'

时间（s）'

title（'

信号波形'

subplot（2,1,2）;

Y1=fft（y1）;

w1=2/n1*（0:

n1-1）;

%设置角频率

plot（w1,abs（Y1））;

%画频谱图

信号频谱'

数字角频率'

幅度'

gridon;

sound（y,fs）;

实验结果：

1、通过观察频谱知，选取音乐信号的频谱集中在0~0.7*pi之间，抽样点数fs=44100;

2、当采样频率问原来0.5（0.5*fs）倍时：

音乐片段音调变得非常低沉，无法辨认原声，播放时间变长；

抽样频率减小，抽样点数不变时，其分辨力增大，记录长度变长，声音失真。

3、当采样频率问原来2（2*fs）倍时：

音乐片段音调变得尖而细，语速变快，播放时间变短；

抽样频率增加，抽样点数不变时，其分辨力下降，记录长度变短，声音失真。

2、音乐信号的抽取（减抽样）

①观察音乐信号频率上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；

②输出减抽样音乐信号的波形和频谱，观察现象，给出理论解释；

③播放减抽样音乐信号，注意抽样率的改变，比较不同抽取间隔下的声音，解释现象。

）

tn1=（0:

figure

plot（tn1,y1）;

原信号波形'

wn1=2/n1*[0:

n1-1];

plot（wn1,abs（Y1））;

原信号频谱'

数字角频率w'

D=2;

%设置抽样间隔

y2=y1（1:

D:

n1）;

%减抽样

n2=length（y2）;

%减抽样后信号长度

t2=（0:

n2-1）/fs;

%设置横坐标

plot（t2,y2）;

%绘制减抽样信号波形图

2:

1减抽样信号波形'

Y2=fft（y2）;

%对y2进行n2点fft谱分析

w2=2/n2*[0:

n2-1];

plot（w2,abs（Y2））;

%绘制减抽样信号频谱图

1减抽样信号频谱'

sound（y2,fs/D）;

实验结果与分析：

1、程序中指标D表示抽样间隔，其值越大，相邻两抽样点之间的距离越远，抽样后漏掉的信息越多，相应的时域信号长度越短；

2、抽样间隔D=1.1时的信号波形及频谱图，抽样频率

大于信号最高频率

的两倍，

满足抽样定理，不会发生混叠。

抽样间隔D越大，抽样率fs越小，抽样后时域信号长度越短

3、抽样间隔D=2时的信号波形及频谱图

抽样间隔D=2的信号波形及频谱图，抽样频率

小于信号最高频率

的两倍，即

<

2

，不满足抽样定理，其频谱图发生混叠，且D越大，混叠越严重，高频成分增加越多，音乐片段音调听起来很沙哑，音调变得很高。

4、抽样间隔D=5的信号波形及频谱图，抽样频率

，不满足抽样定理，其频谱图发生混叠，音乐片段音调听起来很尖锐。

3、音乐信号的AM调制

①观察音乐信号频率上限，选择适当调制频率对信号进行调制（给出高、低两种调制频率）；

②输出调制信号的波形和频谱，观察现象，给出理论解释；

③播放调制音乐信号，注意不同调制频率下的声音，解释现象。

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）;

%信号频谱发生搬移

调制后信号'

w1=2/lc1*[0:

lc1-1];

%设置角频率W

C1=fft（c1）;

subplot（2,1,2）

plot（w1,abs（C1））;

调制后信号的频谱（高频率调制）'

%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（'

y1=a（:

%去单列数据进行分析

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）;

%绘制音乐信号时域波形图

原信号音频'

时间'

subplot（2,2,2）;

plot（w,abs（f1））;

%绘制音乐信号频谱图

频率/pi'

t=[0:

y2=cos（pi*1/2*t）;

%载波函数

y3=y1.*y2'

;

%信号调制

ty3=（0:

length（y3）-1）/fs;

subplot（2,2,3）;

plot（ty3,y3）;

%绘制调制后信号波形图

AM调制音频信号'

f3=fft（y3）;

n2=length（f3）;

subplot（2,2,4）;

plot（w2,abs（f3））;

%绘制调制后信号频谱图

AM调制频谱'

%解调后信号

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）;

解调后的波形'

t'

subplot（3,1,2）;

plot（wy4,abs（Y4））;

解调后波形频谱'

w/pi'

%设计巴特沃斯滤波器进行滤波去噪

[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））;

IIR滤波器频谱'

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）

plot（tm,m）;

IIR滤波器滤波后波形'

plot（wm,abs（M））;

IIR滤波器滤波后频谱'

%矩形窗和布莱克曼窗

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:

M=21184;

fh1=fft（h1,M）;

%矩形窗频谱函数

w=2/M*[0:

M-1];

fh2=fft（h2,M）;

%布莱克曼窗频谱函数

figure（4）

plot（th1,h1）

矩形窗时域'

plot（w,abs（fh1））;

矩形窗频域'

plot（th2,h2）;

布莱克曼窗时域'

plot（w,abs（fh2））;

布莱克曼窗频域'

%滤波处理

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）

plot（ty6,y6）;

矩形窗滤波后音频'

plot（w3,abs（f6））;

矩形窗滤波后频谱'

y7=conv（h2,y4）;

%用布莱克曼窗对调制后的信号进行滤波

f7=fft（y7）;

n5=length（f7）;

ty7=（0:

n5-1）/fs;

w4=2/n5*[0:

n5-1];

%sound（y7,fs）;

plot（ty7,y7）;

布莱克曼窗滤波后音频'

plot（w4,abs（f7））;

布莱克曼窗滤波后频谱'

实验结果与分析：

由调制信号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;

y0=y（:

l=length（y0）;

%加三余弦混合噪声

t0=（0:

l-1）/fs;

d0=[0.05*cos（2*pi*3000*t0）]'

d1=[0.05*cos（2*pi*5000*t1）]'

d2=[0.05*cos（2*pi*8000*t2）]'

noise=d2+d1+d0;

y1=y0+noise;

a=length（noise）;

%绘制三余弦噪声音频图

wa=2/a*[0:

a-1];

Noise=fft（noise）;

plot（noise（1:

150））;

三余弦信号音谱'

%绘制三余弦噪声频谱图

plot（wa,abs（Noise））;

W'

噪声频谱'

w0=2/l*[0:

l-1];

%绘制加噪信号音频

subplot（2,2,4）

plot（w0,abs（Y1））;

加噪信号频谱'

ly1=length（y1）;

ty1=（0:

ly1-1）/fs;

plot（ty1,y1）;

加噪信号音谱'

m=rand（l,1）-0.5;

%产生幅度为0.5的随机信号

y2=m+y0;

%将噪声信号与原声音信号叠加

%sound（y2,fs）;

plot（m（1:

150））

白噪信号音谱'

l=length（y2）;

ty2=（0:

w=2/l*[0:

plot（w,abs（Y2））;

plot（ty2,y2）;

加噪信号幅值'

%设计滤波器进行滤波去噪

[N1,wc1]=buttord（0.04,0.13,1,24）;

%[N1,wc1]=buttord（0.04,0.20,1,18）;

（针对随机白噪声的指标）

figure（3）;

%滤波器频谱

幅度k'

IIR滤波器滤波波形'

m=filter（b,a,y1）;

%用滤波器滤除三余弦噪声

%sound（m,fs）;

%绘制滤波后的波形

t（s）'

信号幅值'

去噪后信号波形'

k=fft（m）;

%滤波后的波形做离散傅里叶变换

w=2*[0:

length（k）-1]/length（k）;

subplot（3,1,3）

plot（w,abs（k））;

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、