使用matlab进行简单音乐合成.docx
《使用matlab进行简单音乐合成.docx》由会员分享，可在线阅读，更多相关《使用matlab进行简单音乐合成.docx（16页珍藏版）》请在冰点文库上搜索。
使用matlab进行简单音乐合成
信号与系统
－综合实验之音乐合成
（1）请根据《东方红》片断的简谱和“十二平均律”计算出该片断中各个乐音的频率，在MATLAB中生成幅度为1、抽样频率为8kHz的正弦信号表示这些乐音。
请用sound函数播放每个乐音，听一听音调是否正确。
最后用这一系列乐音信号拼出《东方红》片断，注意控制每个乐音持续的时间要符合节拍，用sound播放你合成的音乐，听起来感觉如何？
代码如下：
f=8000;
t2=[0:
1/f:
1];
t4=[0:
1/f:
0.5];
t8=[0:
1/f:
0.25];
omg5=523.35;
omg6=587.33;
omg2=392;
omg1=349.23;
omg6l=293.66;
m1=sin（2*pi*omg5*t4）;
m2=sin（2*pi*omg5*t8）;
m3=sin（2*pi*omg6*t8）;
m4=sin（2*pi*omg2*t2）;
m6=sin（2*pi*omg1*t4）;
m7=sin（2*pi*omg1*t8）;
m8=sin（2*pi*omg6l*t8）;
m9=sin（2*pi*omg2*t2）;
m=[m1m2m3m4m6m7m8m9];
sound（m）;
听的时候发现在相邻乐音之间有杂音，这是由于相位不连续造成的。
（2）你一定注意到
（1）的乐曲中相邻乐音之间有“啪”的杂声，这是由于相位不连续产生了高频分量。
这种噪声严重影响合成音乐的质量，丧失真实感。
为了消除它，我们可以用图1.5所示包络修正每个乐音，以保证在乐音的邻接处信号幅度为零。
此外建议用指数衰减的包络来表示。
我采用的是指数衰减的包络。
代码如下：
f=8000;
t2=[0:
1/f:
1];
t4=[0:
1/f:
0.5];
t8=[0:
1/f:
0.25];
omg5=523.35;
omg6=587.33;
omg2=392;
omg1=349.23;
omg6l=293.66;
m1=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m2=exp（-4*t8）.*sin（2*pi*omg5*t8）;
m3=exp（-4*t8）.*sin（2*pi*omg6*t8）;
m4=exp（-1*t2）.*sin（2*pi*omg2*t2）;
m6=exp（-2*t4）.*sin（2*pi*omg1*t4）;
m7=exp（-4*t8）.*sin（2*pi*omg1*t8）;
m8=exp（-4*t8）.*sin（2*pi*omg6l*t8）;
m9=exp（-1*t2）.*sin（2*pi*omg2*t2）;
m=[m1m2m3m4m6m7m8m9];
sound（m）;
第一次我采用的指数衰减没有时间前面的系数，即每个都只乘exp（t）,没有系数；后来根据不同节拍，更改了不同的衰减系数，这样声音听起来感觉更加圆润。
（3）请用最简单的方法将
（2）中的音乐分别升高和降低一个八度。
（提示：
音乐播放的时间可以变化）再难一些，请用resample函数（也可以用interp和decimate函数）将上述音乐升高半个音阶。
（提示：
视计算复杂度，不必特别精确）
答：
最简单的方法是直接更改抽样频率f。
将f从8K改为4K，则升高一个八度，并且播放速度增快了一倍；将f从8k改为16k，则降低一个八度，速度也变慢了一倍。
升高半个音阶，只须在
（2）代码最后加一句resample（m,1000,1059）即可。
f=8000;%改为4000或者16000
t2=[0:
1/f:
1];
t4=[0:
1/f:
0.5];
t8=[0:
1/f:
0.25];
omg5=523.35;
omg6=587.33;
omg2=392;
omg1=349.23;
omg6l=293.66;
m1=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m2=exp（-4*t8）.*sin（2*pi*omg5*t8）;
m3=exp（-4*t8）.*sin（2*pi*omg6*t8）;
m4=exp（-1*t2）.*sin（2*pi*omg2*t2）;
m6=exp（-2*t4）.*sin（2*pi*omg1*t4）;
m7=exp（-4*t8）.*sin（2*pi*omg1*t8）;
m8=exp（-4*t8）.*sin（2*pi*omg6l*t8）;
m9=exp（-1*t2）.*sin（2*pi*omg2*t2）;
m=[m1m2m3m4m6m7m8m9];
resample（m,1000,1059）;
sound（m）;
（4）试着在
（2）的音乐中增加一些谐波分量，听一听音乐是否更有“厚度”了？
注意谐波分量的能量要小，否则掩盖住基音反而听不清音调了。
（如果选择基波幅度为1，二次谐波幅度0:
2，三次谐波幅度0:
3，听起来像不像象风琴？
）
代码如下：
f=8000;
t2=[0:
1/f:
1];
t4=[0:
1/f:
0.5];
t8=[0:
1/f:
0.25];
omg5=523.35;
omg6=587.33;
omg2=392;
omg1=349.23;
omg6l=293.66;
m1=exp（-2*t4）.*sin（2*pi*omg5*t4）+0.2*exp（-2*t4）.*sin（2*pi*2*omg5*t4）+0.3*exp（-2*t4）.*sin（2*pi*3*omg5*t4）;
m2=exp（-4*t8）.*sin（2*pi*omg5*t8）+0.2*exp（-4*t8）.*sin（2*pi*2*omg5*t8）+0.3*exp（-4*t8）.*sin（2*pi*3*omg5*t8）;
m3=exp（-4*t8）.*sin（2*pi*omg6*t8）+0.2*exp（-4*t8）.*sin（2*pi*2*omg6*t8）+0.3*exp（-4*t8）.*sin（2*pi*3*omg6*t8）;
m4=exp（-1*t2）.*sin（2*pi*omg2*t2）+0.2*exp（-1*t2）.*sin（2*pi*2*omg2*t2）+0.3*exp（-1*t2）.*sin（2*pi*3*omg2*t2）;
m6=exp（-2*t4）.*sin（2*pi*omg1*t4）+0.2*exp（-2*t4）.*sin（2*pi*2*omg1*t4）+0.3*exp（-2*t4）.*sin（2*pi*3*omg1*t4）;
m7=exp（-4*t8）.*sin（2*pi*omg1*t8）+0.2*exp（-4*t8）.*sin（2*pi*2*omg1*t8）+0.3*exp（-4*t8）.*sin（2*pi*3*omg1*t8）;
m8=exp（-4*t8）.*sin（2*pi*omg6l*t8）+0.2*exp（-4*t8）.*sin（2*pi*2*omg6l*t8）+0.3*exp（-4*t8）.*sin（2*pi*3*omg6l*t8）;
m9=exp（-1*t2）.*sin（2*pi*omg2*t2）+0.2*exp（-1*t2）.*sin（2*pi*2*omg2*t2）+0.3*exp（-1*t2）.*sin（2*pi*3*omg2*t2）;
m=[m1m2m3m4m6m7m8m9];
sound（m）;
加入谐波分量后，音色有所变化，感觉更加清脆一些。
（5）自选其它音乐合成，例如贝多芬第五交响乐的开头两小节。
我选取的是《晴天》的第一句
代码如下：
f=8000;
t2=[0:
1/f:
1];
t4=[0:
1/f:
0.5];
t8=[0:
1/f:
0.25];
t=[0:
1/f:
0.125];
omg1=392;
omg2=440;
omg3=493.88;
omg4=523.25;
omg5=587.33;
omg6=659.25;
omg7=698.45;
omg5l=293.66;
m0=0;
m1=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m2=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m3=exp（-2*t4）.*sin（2*pi*omg1*t4）;
m4=exp（-1*t2）.*sin（2*pi*omg1*t2）;
m5=exp（-2*t4）.*sin（2*pi*omg2*t4）;
m6=exp（-2*t2）.*sin（2*pi*omg3*t2）;
m7=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m8=exp（-2*t4）.*sin（2*pi*omg5*t4）;
m9=exp（-2*t4）.*sin（2*pi*omg1*t4）;
m10=exp（-2*t4）.*sin（2*pi*omg1*t4）;
m11=exp（-4*t8）.*sin（2*pi*omg2*t8）;
m12=exp（-4*t8）.*sin（2*pi*omg3*t8）;
m13=exp（-4*t8）.*sin（2*pi*omg2*t8）;
m14=exp（-4*t8）.*sin（2*pi*omg1*t8）;
m15=exp（-2*t2）.*sin（2*pi*omg5l*t2）;
m=[m0m1m2m3m4m5m6m7m8m9m10m11m12m13m14m15];
sound（m）;
（6）先用wavread函数载入光盘中的fmt.wav文件，播放出来听听效果如何？
是否比刚才的合成音乐真实多了？
x=wavread（'fmt.wav'）;
sound（x）;
（7）你知道待处理的wave2proc是如何从真实值realwave中得到的么？
这个预处理过程可以去除真实乐曲中的非线性谐波和噪声，对于正确分析音调是非常重要的。
提示：
从时域做，可以继续使用resample函数。
realwave中的波形有十个周期，要除去其中的噪声可以采用时域求均值的方法。
步骤如下：
首先用resample函数将其采样率增大为十倍，再等分十份取平均，然后重复这个平均后的波形十次，还原其长度，最后在用resample函数还原采样率到1/10。
代码如下：
load（'Guitar.MAT'）;
wavetemp=zeros（length（realwave）,1）;
waveresampled=resample（realwave,10,1）;
forn=1:
10
wavetemp=wavetemp+waveresampled（（n-1）*length（realwave）+1:
n*length（realwave））/10;
end
mywave2proc=repmat（wavetemp,10,1）;
mywave2proc=resample（mywave2proc,1,10）;
figure;
subplot（3,1,1）;plot（realwave）;
subplot（3,1,2）;plot（wave2proc）;
subplot（3,1,3）;plot（mywave2proc）;
（8）这段音乐的基频是多少？
是哪个音调？
请用傅里叶级数或者变换的方法分析它的谐波分量分别是什么。
提示：
简单的方法是近似取出一个周期求傅里叶级数但这样明显不准确，因为你应该已经发现基音周期不是整数（这里不允许使用resample函数）。
复杂些的方法是对整个信号求傅里叶变换（回忆周期性信号的傅里叶变换），但你可能发现无论你如何提高频域的分辨率，也得不到精确的包络（应该近似于冲激函数而不是sinc函数），可选的方法是增加时域的数据量，即再把时域信号重复若干次，看看这样是否效果好多了？
请解释之。
答：
共计221-2=219个周期；基因频率约为8000/（219/9）=328.7671Hz.
用FFT对该信号做DFT变换：
load（'Guitar.MAT'）;
my_wave2proc=repmat（wave2proc,25,1）;
Fs=8000;
Length=length（my_wave2proc）;
NFFT=2^nextpow2（Length）;
Y=fft（my_wave2proc,NFFT）/Length;
amplitude=2*abs（Y（1:
NFFT/2+1））;
frequency=Fs/2*linspace（0,1,NFFT/2+1）;
plot（frequency,amplitude）;
[max_rate_of_grade,max_position]=max（amplitude（1:
100））;
frequency（max_position）
算得基频为329.1016Hz,和前面通过周期估算得到的很相近。
对照可得，为C大调mi。
（9）再次载入fmt.wav，现在要求你写一段程序，自动分析出这段乐曲的音调和节拍！
如果你觉得太难就允许手工标定出每个音调的起止时间，再不行你就把每个音调的数据都单独保存成一个文件，然后让MATLAB对这些文件进行批处理。
注意：
不允许逐一地手工分析音调。
编辑音乐文件，推荐使用\CoolEdit"编辑软件。
此问参考了学长的版本。
看了版本才明白了算法。
主要分2个步骤：
1、分割出单个音符，计算每个音符的时间，2、计算每个音符的频率，把频率转化为音名。
分割音符
即将一个一个的音符片段从整个曲子中切下。
音符的起始都伴随时域上瞬时能量的激增，即时域上波形幅度突然变大。
利用这一点特征对音符进行提取。
选择小步幅对fmt.wav进行扫描，步幅选为0.01秒最为适宜。
具体方法如下：
在0.01秒的范围内找极大点，而且是后续能量连续下降的极大点，这个极大点如果满足以下三个条件，则为有效的音符的起始点：
1、这个极大点比前一段的极大点能量高出80%；Pw>1.8*Pwp
2、这个极大点比整个乐曲的平均能量高；St>Stp+1000
3、距前面最近的一个音符的起始点时间差1/8秒以上。
Pw>Avg
分割的结果：
一共分割得30个音，和CoolEdit的分析基本一致。
过度部分
每个音符的持续时间Time=diff（[StartPointLen]）/8000;
用于傅立叶变换的区间EndPoint=[StartPoint（2:
Num）Len]-200;
能量修正（归一）Power=Power/max（Power）;
用于傅立叶变换的区间长度STime=Time-200/8000;
计算频率并转化
对每个音符内部的一个区间（我选这个音符的起始到结束前的0.025秒，即StartPoint到EndPoint）做快速傅立叶变换，然后在频谱中选出基频。
关于自动寻找基频：
一般的，基波不是幅度最大的频率，因为有很多音调的谐波，甚至非线性谐波的幅度大于基波。
一般来说，一次谐波的幅度不太大，不会大于基波幅度的2倍，这样，在有可能出现基波的地方，利用每个点的幅度值除以该点的频率，也就是每个点对于原点的斜率做为判定基波的标准，经过这样的简化处理，找到斜率最大的点，基本上就能够判定基波出现在什么位置。
FO=F_fft./omg;%每个点的幅度值除以该点的频率
[M1M2]=max（FO）;%找到斜率最大的点，判定为基波
Freq（p）=omg（M2）/2/pi;%记下频率值
Name（p）=Freq2Name（Freq（p））;%把频率值转化为音名
（关于函数Freq2Name（Freq）：
把频率与十几个标准音的频率作比较，选出最接近的一个，作为这个频率对应的音名。
）
display（'Score：
'）;display（Name）;display（Time）;%输出结果：
#
音名
时长
#
音名
时长
#
音名
时长
#
音名
时长
1
-6
0.1676
9
-4
0.5044
17
-6
0.3890
25
1
0.4718
2
-6
1.4866
10
-4
0.7484
18
5
0.2687
26
-4
0.5336
3
-7
0.4764
11
-7
1.2373
19
4
0.2305
27
-6
0.4627
4
-6
0.4591
12
3
0.7869
20
3
0.2010
28
-7
0.5771
5
2
0.4584
13
-6
0.2236
21
2
0.2610
29
-6
0.6410
6
3
0.4500
14
-6
0.4713
22
1
0.4731
30
-5
1.3974
7
-5
0.4660
15
-6
0.4728
23
-7
0.4437
8
-4
0.4469
16
6
0.4979
24
2
0.5484
如此，便完成了乐曲的音调和节拍的自动分析。
把此结果用函数CreateMusic处理并回放，听到的音乐是很接近原来的fmt.wav的，仅仅有两个左右的音符有偏差。
代码如下：
fmt=wavread（'fmt.wav'）;
%sound（fmt,8000）
Len=length（fmt）;
Avg=norm（fmt）/sqrt（Len）;
Num=0;Step=80;Pwp=0;Stp=-1000;St=1;
forn=1:
Step:
Len-Step
M=0;St=n;
whileM~=1&St+Step-1<=Len
[PwM]=max（fmt（St:
St+Step-1））;
St=St+M-1;
end
ifPw>1.8*Pwp&St>Stp+1000&Pw>Avg
Num=Num+1;
StartPoint（Num）=St;
Power（Num）=Pw;
Stp=St;
end
Pwp=Pw;
end
Time=diff（[StartPointLen]）/8000;
STime=Time-200/8000;
EndPoint=[StartPoint（2:
Num）Len]-200;
Power=Power/max（Power）;
figure;holdon;plot（fmt）;
axis（[0length（fmt）00.6]）;
title（'Òô·û·Ö¸î'）;
xlabel（'Time'）;ylabel（'Signal'）;
forp=1:
Num
plot（[StartPoint（p）StartPoint（p）],[00.6],'r*-'）;
end
FreqLimit=170;
forp=1:
Num
f=fmt（StartPoint（p）:
EndPoint（p）-1）;
[F_fft,t,omg]=FastFourier（f）;
omg（floor（STime（p）*（4000-FreqLimit））:
ceil（STime（p）*（4000+FreqLimit）））=1e5;
FO=F_fft./omg;
[M1M2]=max（FO）;
Freq（p）=omg（M2）/2/pi;
Name（p）=Freq2Name（Freq（p））;
end
display（'Score£º'）;display（Name）;display（Time）;
UnknownMusic=[Name;Time;Power];
Music=CreateMusic（'C',UnknownMusic,'Exp','Guitar+'）;
PlayAndPlot（Music）;
（10）用（7）计算出来的傅里叶级数再次完成第（4）题，听一听是否像演奏fmt.wav的吉他演奏出来的？
答：
在（7）中计算出来的傅里叶级数中取前四项，后面的谐波并不显著将其略去。
基波和2、3、4次谐波的系数为如下，将其归一后再合成即可。
合成的方法参考了学长的报告。
music_score=[55621162;
000000-10;
10.50.5210.50.52];
time=5;
harmonic_coefficient=[0.05380.07530.04640.0482];
harmonic_coefficient=harmonic_coefficient/norm（harmonic_coefficient）;
diationic_scale=[220*2.^（（[3578101214]+5）/12）];
t=[0:
1/8000:
time];
music_wave=zeros（1,length（t））;
time_start=0;
speed=0.5;
forn=1:
length（music_score）
frequency=diationic_scale（music_score（1,n））*2^music_score（2,n）;
time_end=time_start+music_score（3,n）*speed;
envelope=GetEnvelope3（time,time_start,time_end）;
form=1:
4
music_wave=music_wave+harmonic_coefficient（m）*sin（2*pi*m*frequency*t）.*（t>=time_start&tend
time_start=time_end;
end
music_wave=music_wave/max（music_wave）;
sound（music_wave,8000）;
（11）也许（9）还不是很像，因为对于一把泛音丰富的吉他而言，不可能每个音调对应的泛音数量和幅度都相同。
但是通过完成第（8）题，你已经提取出fmt.wav中的很多音调，或者说，掌握了每个音调对应的傅里叶级数，大致了解了这把吉他的特征。
现在就来演奏一曲《东方红》吧。
提示：
如果还是音调信息不够，那就利用相邻音调的信息近似好了，毕竟可以假设吉他的频响是连续变化的。
答：
和上题类似。
music_score=[55621162;
000000-10;
10.50.5210.50.52];
time=5;
load（'HarmonicCoefficients.mat'）;
diationic_scale=[220*2.^（（[3578101214]+5）/12）];
t=[0:
1/8000:
升级会员 