Matlab音乐合成实验报告.docx
《Matlab音乐合成实验报告.docx》由会员分享，可在线阅读，更多相关《Matlab音乐合成实验报告.docx（25页珍藏版）》请在冰点文库上搜索。
Matlab音乐合成实验报告
音乐合成实验
摘要：
本文共有三大部分：
第一部分，简单的音乐合成；第二部分，用傅里叶变换分析音乐；第三部分，基于傅里叶级数的音乐合成。
由潜入深，一步一步分析了用MATLAB进行音乐合成的过程。
通过本实验达到了加深对傅里叶级数和傅里叶分析的理解，熟悉对MATLAB基本使用的目标。
第一部分简单的合成音乐
1.1合成《东方红》
根据《东方红》第一小节的简谱和十二平均律计算出该小节每个乐音的频率，在MATLAB中生成幅度为1，抽样频率为8kHz的正弦信号表示这些乐音，用sound播放合成的音乐
由图可知《东方红》的曲调定为F，即1=F，对应的频率为349.23Hz，据此可以计算出其他乐音的频率，例如5对应的频率为
，一次类推计算出第一小节各乐音对应的频率为：
乐音
5
5
6
2
1
1
6
2
频率
523.25
523.25
587.33
392
349.23
349.23
293.66
392
在确定了各乐音的频率之后需要确定每个乐音的持续时间。
每小节有两拍，一拍的时间是0.5s，因此各乐音的持续时间为：
乐音
5
5
6
2
1
1
6
2
时间
0.5
0.25
0.25
1
0.5
0.25
0.25
1
而在MATLAB中表示乐音所用的抽样频率为fs=8000Hz，也就是所1s钟内有8000个点，抽样点数的多少就可表示出每个乐音的持续时间的长短。
用一个行向量来存储这段音乐对应的抽样点，在用sound函数播放即可。
根据以上分析在MATLAB中编写如下程序：
sound_1_1.m
clear;clc;
fs=8000;%抽样频率
f=[523.25523.25587.33392349.23349.23293.66392];
%各个乐音对应的频率
time=fs*[1/2,1/4,1/4,1,1/2,1/4,1/4,1];%各个乐音的抽样点数
N=length（time）;%这段音乐的总抽样点数
east=zeros（1,N）;%用east向量来储存抽样点
n=1;
fornum=1:
N%利用循环产生抽样数据，num表示乐音编号
t=1/fs:
1/fs:
time（num）/fs;%产生第num个乐音的抽样点
east（n:
n+time（num）-1）=sin（2*pi*f（num）*t）;
%抽样点对应的幅值
n=n+time（num）;
end
sound（east,8000）;%播放音乐
在MATLAB中运行sound_1_1.m，播放出了《东方红》的第一段，但是可以听出效果很不好，只能听出具有《东方红》的调子而已。
1.2除噪音，加包络
在1.1中听到有“啪”的杂声，下面通过加包络来消噪音。
最简单的包络为指数衰减。
最简单的指数衰减是对每个音乘以
因子，在实验中首先加的是
的衰减，这种衰减方法使用的是相同速度的衰减，但是发现噪音并没有完全消除，播放的音乐效果不是很好，感觉音乐起伏性不强。
于是采用不同速度的衰减，根据乐音持续时间的长短来确定衰减的快慢，乐音持续时间越长，衰减的越慢，持续时间越短，衰减的越快。
在1.1程序的基础上加上包络，编写如下程序：
sound_1_21.m
clear;clc;
fs=8000;%抽样频率
f=[523.25523.25587.33392349.23349.23293.66392];
%各个乐音对应的频率
time=fs*[1/2,1/4,1/4,1,1/2,1/4,1/4,1];%各个乐音的抽样点数
N=length（time）;%这段音乐的总抽样点数
east=zeros（1,N）;%用east向量来储存抽样点
n=1;
fornum=1:
N%利用循环产生抽样数据，num表示乐音编号
t=1/fs:
1/fs:
time（num）/fs;%产生第num个乐音的抽样点
G=zeros（1,time（num））;%G为存储包络数据的向量
G（1:
time（num））=exp（1:
（-1/time（num））:
1/8000）;
%产生包络点
east（n:
n+time（num）-1）=sin（2*pi*f（num）*t）.*G（1:
time（num））;
%给第num个乐音加上包络
n=n+time（num）;
end
sound（east,8000）;%播放
plot（east）;
播放后可以听出噪音已经消除，同时因为不同时长的乐音衰减的快慢不一样，音乐听起来更有起伏感，下图是加包络后的east图像。
更科学的包络如下图所示，每个乐音都经过冲激、衰减、持续、消失四个阶段。
由上图可以看出这个包络是四段直线段构成的，因此只要确定了每段线段的端点，即可用端点数据写出直线方程，因为直线方程可以用通式写出（我用的是斜截式），因此这段包络可以用简单的循环来完成。
例如认为包络线上的数据如下图所示：
据此在MATLAB中编写如下程序：
sound_1_22.m
clear;clc;
fs=8000;%抽样频率
f=[523.25523.25587.33392349.23349.23293.66392];
%各个乐音对应的频率
time=fs*[1/2,1/4,1/4,1,1/2,1/4,1/4,1];%各个乐音的抽样点数
N=length（time）;%这段音乐的总抽样点数
east=zeros（1,N）;%用east向量来储存抽样点
n=1;
fornum=1:
N%利用循环产生抽样数据，num表示乐音编号
t=1/fs:
1/fs:
（time（num））/fs;%产生第num个乐音的抽样点
P=zeros（1,time（num））;%P为存储包络数据的向量
L=（time（num））*[01/5333/1000333/5001];
%包络线端点对应的横坐标
T=[01.5110];%包络线端点对应的纵坐标
s=1;
b=1:
1:
time（num）;%产生包络线抽样点
fork=1:
4P（s:
L（k+1）-1）=（T（k+1）-T（k））/（L（k+1）-L（k））*（b（s:
L（k+1）-1）-L（k+1）*ones（1,L（k+1）-s））+T（k+1）*ones（1,L（k+1）-s）;
%包络线直线方程通式
s=L（k+1）;
endeast（n:
n+time（num）-1）=sin（2*pi*f（num）*t）.*P（1:
time（num））;
%给第num个乐音加上包络
n=n+time（num）;
end
sound（east,8000）;
plot（east）;
运行得到的图像为：
下图是两个乐音交接处的局部放大图，可以看到前一个乐音一直衰减到0，后一个乐音从0开始增加，因此消除了噪音。
若不需要每个音都衰减到0，例如只需衰减到持续阶段幅值的20%，对程序做简单的修改即可，将T=[01.5110]改为T=[0.21.5110.2]运行得到的结果为：
由图可见，每个乐音都是衰减到一较小值而不是0，也能消除噪音，同时音乐听起来更加连续。
1.3改变程序，实现1.2中的音乐升高和降低一个八度
升高一个八度即每个乐音的频率都提高一倍，变为原来的2被；降低一个八度即每个乐音的频率都减小一倍，变为原来的1/2。
因此最简单的办法是将存储乐音频率的向量每个元素改变为2或1/2倍。
即将程序中的f=[523.25523.25587.33392349.23349.23293.66392];改为
f=[523.25523.25587.33392349.23349.23293.66392]*2;或
f=[523.25523.25587.33392349.23349.23293.66392]/2;
将上述音乐上高半个音阶，即将频率变为原来的
（1.06）倍，可以利用resamlpe函数对原来的数据点进行重采样来实现
east=resample（east,100,106）;
因为resample进行重新采样后会使每个乐音的持续时间改变，但是因为升高半个音阶，频率改变不大，所以每个音的持续时间是基本不变的。
1.4在1.2的音乐中加入谐波
在1.2的音乐中加上二、三、四次谐波，基波幅度为1，高次谐波幅度分别为0.2、0.3、0.1。
只需将1.2程序改为
sound_1_4.m
clear;clc;
fs=8000;%抽样频率
f=[523.25523.25587.33392349.23349.23293.66392];
%各个乐音对应的频率
time=fs*[1/2,1/4,1/4,1,1/2,1/4,1/4,1];%各个乐音的抽样点数
N=length（time）;%这段音乐的总抽样点数
east=zeros（1,N）;%用east向量来储存抽样点
n=1;
fornum=1:
N%利用循环产生抽样数据，num表示乐音编号
t=1/fs:
1/fs:
（time（num））/fs;%产生第num个乐音的抽样点
P=zeros（1,time（num））;%P为存储包络数据的向量
L=（time（num））*[01/5333/1000333/5001];
%包络线端点对应的横坐标
T=[01.5110];%包络线端点对应的纵坐标
s=1;
b=1:
1:
time（num）;%产生包络线抽样点
fork=1:
4P（s:
L（k+1）-1）=（T（k+1）-T（k））/（L（k+1）-L（k））*（b（s:
L（k+1）-1）-L（k+1）*ones（1,L（k+1）-s））+T（k+1）*ones（1,L（k+1）-s）;
%包络线直线方程通式
s=L（k+1）;
end
m=[10.30.2];%波形幅值矩阵
ss=zeros（1,length（t））;
fori=1:
length（m）
ss=ss+m（i）*sin（2*i*pi*f（num）*t）;%加谐波
end
east（n:
n+time（num）-1）=ss.*P（1:
time（num））;
%给第num个乐音加上包络
n=n+time（num）;
end
sound（2*east,8000）;
plot（east）;
即可，加颜色部分为修改的部分，加上谐波后音乐效果变得更好了。
1.5自选音乐合成——《两只老虎》
两只老虎
曲调为C，因此可以得到每个乐音对应的频率分别为：
1
2
3
4
5
262.63
293.66
329.63
349.23
392
每小节有四拍，一拍的时间是0.5s，因此各乐音的持续时间为：
乐音
1
2
3
1
1
2
3
1
时间
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
乐音
3
4
5
3
4
5
时间
0.5
0.5
1
0.5
0.5
1
0.25
sound_1_5.m
clear;clc;
fs=8000;%抽样频率
f=[262.63293.66329.63262.63262.63293.66329.63262.63329.63349.23392];
%各个乐音对应的频率
time=fs*[1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2,1,1/2,1/2,1];
%各个乐音的抽样点数
N=length（time）;%这段音乐的总抽样点数
east=zeros（1,N）;%用east向量来储存抽样点
n=1;
fornum=1:
N%利用循环产生抽样数据，num表示乐音编号
t=1/fs:
1/fs:
time（num）/fs;%产生第num个乐音的抽样点
G=zeros（1,time（num））;%G为存储包络数据的向量
G（1:
time（num））=exp（1:
（-1/time（num））:
1/8000）;
%产生包络点
east（n:
n+time（num）-1）=sin（2*pi*f（num）*t）.*G（1:
time（num））;
%给第num个乐音加上包络
n=n+time（num）;
end
sound（east,8000）;%播放
plot（east）;
第二部分用傅里叶变换分析音乐
2.1载入fmt.wav并播放
利用wavread函数载入，用sound函数播放，程序如下:
sound_2_1.m
wave=wavread（'fmt.wav'）;
sound（wave）
这段音乐听起来比之前合成的音乐更加真实，因为里边含有丰富的谐波。
2.2载入文件Guitar.mat，处理原始数据realwave
载入文件Guitar.mat，分析wave2proc是怎么由realwave得到的。
利用
loadGuitar.mat;载入并用plot函数将realwave、wave2proc分别画出，得到以下两幅图
可以看到，wave2proc比realwave的周期性好得多，去掉了非线性谐波和噪声。
在时域做，从图上可以看到，realwave的数据大约是10个周期的共243个数据，因此可以用resample函数对realwave进行重新采样，将采样点提高到250个，那么重采样后每个周期有25个点，将这25个点对应相加求平均值后得到一个周期的值，因为进行了平均，减小了非线性谐波和噪音，然后将这25个数据延托成十个周期即250个点，在利用resample函数对得到的函数重新采样将采样点数恢复到243个。
根据以上分析，编写实现这个思路的程序如下：
sound_2_2.m
clear;clc;
loadGuitar.mat;
wave=resample（realwave,250,243）;%重采样，将点数变为250
w=zeros（1,25）;
fori=1:
25
fork=0:
9
w（i）=w（i）+wave（25*k+i）;%10个周期的对应点分别求和
end
end
w=w/10;%取平均值
wave2=repmat（w,1,10）;%将1个周期的10个点延拓至250个点
wave2=resample（wave2,243,250）;%重采样，将点数变回243
holdon,plot（wave2,'r'）,holdoff;%将处理后的数据绘出，红色
holdon,plot（wave2proc）;%将所给的数据绘出，蓝色
运行后的结果为：
由图可见，两组数据重合的很好，说明这种方法是很不错的方法。
2.3分析wave2proc的基波和谐波
为了分析wave2proc的基波和谐波，可以对wave2proc进行傅里叶变换，得到wave2proc的幅值谱，在频谱图上的第一个突出的波峰对应的频率即为wave2proc基频，利用helpfft学习了MATLAB中快速傅里叶变换函数fft的用法，编写了如下程序：
clear;clc;
loadGuitar.mat;
fs=8000;
NFFT=2^nextpow2（length（wave2proc））;
Y=fft（wave2proc,NFFT）/length（wave2proc）;
g=fs/2*linspace（0,1,NFFT/2+1）;
plot（g,2*abs（Y（1:
NFFT/2+1）））
运行后得到的结果为
虽然从图上可以大概看出包络，但是非常不明显，假如提高频域的抽样频率，例如将抽样频率由NFFT=2^nextpow2（length（wave2proc））改为
NFFT=8^nextpow2（length（wave2proc））得到的结果如下;
由图可见虽然频域的抽样频率提高了很多，但是得到的包络依然不精确，这是因为wave2proc是周期函数，但是现在的wave2proc只有243个数据点，并不能非常明显的体现出其周期性，因此它的幅值谱的离散化程度不高，虽然提高了频域的抽样频率，但是wave2proc数据点的周期性并没有增加，所以要显示出离散化程度高的幅值谱，就要增加wave2proc的周期性，即让wave2proc在时域重复多次后在进行傅里叶变换。
利用repmat函数可以将wave2proc在时域重复。
将程序修改为
sound_2_3.m
clear;clc;
loadGuitar.mat;
fs=8000;
wave2proc=repmat（wave2proc,20,1）;%将wave2proc重复20次
NFFT=2^nextpow2（length（wave2proc））;
Y=fft（wave2proc,NFFT）/length（wave2proc）;
g=fs/2*linspace（0,1,NFFT/2+1）;
plot（g,2*abs（Y（1:
NFFT/2+1）））
运行后得到的幅值谱为
可以看出幅值谱的离散化程度已经非常高了。
由图读出wave2proc的基频为329.1Hz，幅值为0.05401，高次谐波幅值分别为：
谐波
2
3
4
5
6
7
8
9
幅值
0.07676
0.04841
0.0519
0
0.005709
0.01923
0.006741
0.007326
2.4自动分析fmt.wav的音调和节拍
思路分析：
将fmt.wav导入后得到的是一个向量，它包含了这段音乐的所有信息，要自动分析这段音乐的音调就需要将每个音调对应的点进行傅里叶变换得到其幅值谱，在幅值谱上找到第一个幅值较大的极大值点，该点对应的就是该音调的基频，得到基频后就可以得到高次谐波的幅值。
为了使对每个音调进行傅里叶变换后得到的幅值谱离散程度高，应该将每个音调的数据在时域上重复多次，由于这些点都是直接采集的为做处理的点，因此其重复次数应该足够大才能体现出较强周期性，本实验采用重复1000次，虽然重复次数越多越好，但是次数太大，程序运行的速度会大大降低。
这里边还有两个关键点：
第一，在从幅值谱上找基频时，因为图上的极大值点很多，怎么能让程序自动确定出准确的基频。
第二，在程序找到了基频之后，再由基频去获取高次谐波的幅值时需要有一定的容错能力，例如若基频为200Hz，幅值为1，那么对应的二次谐波的频率为400Hz，但是很可能恰好幅值谱上400Hz处的幅值为0.01，但是401Hz处的幅值为0.2，这时实际上的二次谐波应该为401Hz，但是若没有给基频一个容错范围，显然找到的二次谐波的幅值是不正确的。
针对以上提出的两个关键点，我找到了两条有针对性的解决办法。
对于第一点，因为幅值谱上极大值点的幅值足够大才能将其定位基频，因此在分析了几个音调后发现基频处的幅值都在0.025以上，因此将基频处的限定条件改为幅值大于0.025的，但是在运行后发现，有几个音调没有分析出来，说明它们的基频幅值小于0.025，其实可以观察一下fmt.wav的波形就会发现，有几段的整体幅值很小，因此基频幅值小，于是又在加上限定条件，若所有点的幅值都小于0.02，那么再用0.015作为幅值的限定条件继续找，这样就将剩下的音调基频也确定出来了。
对于上述的第二点，可以将确定出的基频的误差设为+-1Hz，例如程序确定的基频为200Hz，实际的基频应该在（200-1）到（200+1）之间，那么k次谐波对应的频率范围是k*（200-1）到k*（200+1），在这个区间中继续找幅值的极大值点就是k次谐波对应点。
根据以上思路，下面开始编写用于分析一个音调频率的函数analysis。
在取谐波幅值时，幅值小于基波幅值5%的谐波认为其幅值为0，最终谐波的幅值用归一化后的数据表示。
每一步的详细思路见注释。
analysis.m
function[y1y2]=analysis（w,a）
%设有两个返回值，y1返回频率，y2返回幅值，两个变量，w为待分析数组，a为数据的抽样频率
fs=a;%傅里叶变换的抽样频率
y1=zeros（1,7）;%求最大7次谐波，因此定义1*7矩阵
y2=zeros（1,7）;
NFFT=2^nextpow2（length（w））;
Y=fft（w,NFFT）/length（w）;
g=fs/2*linspace（0,1,NFFT/2+1）;
p=2*abs（Y（1:
NFFT/2+1））;
plot（g,p）%以上为傅里叶变换部分
d=floor（NFFT/fs）;%将误差1Hz化成对应的点数
fork=2:
length（p）-1
if（p（k）>0.02）&&（p（k）>p（k-1））&&（p（k）>p（k+1））%寻找基频的条件
y1
（1）=g（k）;%存储基频
y2
（1）=p（k）;%存储基波幅值
break
elseif（p（k）>0.015）&&（p（k）>p（k-1））&&（p（k）>p（k+1））
%若未找到基频，将幅值限制改为0.015，继续寻找
y1
（1）=g（k）;
y2
（1）=p（k）;
end
end
fort=2:
7
fori=t*（k-d）:
t*（k+d）%在误差允许的范围内寻找t次谐波点
if（p（k）>0.02）&&（p（i）>0.05*p（k））&&（p（i）>p（i-1））&&（p（i）>p（i+1））
y2（t）=p（i）/y2
（1）;%谐波幅值归一化
y1（t）=g（i）;
break
elseif（p（k）>0.015）&&（p（i）>0.05*p（k））&&（p（i）>p（i-1））&&（p（i）>p（i+1））
y2（t）=p（i）/y2
（1）;
y1（t）=g（i）;%这一段循环是与找基频条件相对应的获取高次谐波幅值
end
end
end
在编写完分析函数analysis后即可编写自动分析的主程序。
首先在CooolEdit中手动标定音调交界处的时间节点，得到了time向量
time=floor（[0.0960.2671.7672.2342.7063.1463.606...
4.0564.5205.0305.7495.9787.0157.7097.923...
8.0288.4908.9599.4549.85210.12510.356...
10.56510.82211.29211.74112.28412.741...
13.26913.758
升级会员 