Matlab音乐合成实验报告

1、音乐合成实验目录音乐合成实验 摘要：第一部分简单的合成音乐1.1 合成东方红1.2 除噪音，加包络1.3改变程序，实现1.2中的音乐升高和降低一个八度1.4在1.2的音乐中加入谐波1.5自选音乐合成两只老虎 第二部分用傅里叶变换分析音乐2.1载入 fmt.wav 并播放 2.2载入文件Guitar.mat，处理原始数据realwave2.3分析wave2proc的基波和谐波2.4自动分析fmt.wav的音调和节拍第三部分基于傅里叶级数的音乐合成 3.1用2.3分析出来的结果重新加谐波 3.2通过2.4提取的吉他音调信息弹奏东方红.实验收获摘要：本文共有三大部分：第一部分，简单的音乐合成；第二部

2、分，用傅里叶变换分析音乐； 第三部分，基于傅里叶级数的音乐合成。由潜入深，一步一步分析了用MATLABS行音乐合成的过程。通过本实验达到了加深对傅里叶级数和傅里叶分析的理解，熟悉对MATLABS本使用的目标第一部分简单的合成音乐1.1 合成东方红根据东方红第一小节的简谱和十二平均律计算出该小节每个乐音的频率，在MATLAB中生成幅度为1，抽样频率为8kHz的正弦信号表示这些乐音，用sou nd播放合成的音乐由图可知东方红的曲调定为 F,即 仁F,对应的频率为349.23Hz,据此可以计算出7/12其他乐音的频率，例如5对应的频率为f5 = 349.23 2= 523.25，次类推计算出第一小节

3、各乐音对应的频率为：乐音 55621162频率523.25523.25587.33392349.23349.23293.66392在确定了各乐音的频率之后需要确定每个乐音的持续时间。每小节有两拍，一拍的时间是0.5s，因此各乐音的持续时间为：乐音55621162时间0.50.250.2510.50.250.251而在MATLA中表示乐音所用的抽样频率为fs=8000Hz,也就是所1s钟内有8000个点，抽 样点数的多少就可表示出每个乐音的持续时间的长短。 用一个行向量来存储这段音乐对应的 抽样点，在用sou nd函数播放即可。根据以上分析在MATLA中编写如下程序：soun d_1_1.mcl

4、ear;clc;fs=8000;%由样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率 time=fs*1/2,1/4,1/4,1,1/2,1/4,1/4,1;%各个乐音的抽样点数N=le ngth(time); %这段音乐的总抽样点数 east=zeros(1,N); %用east向量来储存抽样点n=1;for num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/fs:1/fs:time(num)/fs;%产生第nun个乐音的抽样点east( n:n+time( nu m)-1)=s in (2*pi

5、*f( nu m)*t);%由样点对应的幅值n=n+time( nu m);endsou nd(east,8000);%播放音乐在MATLA中运行sound_1_1.m,播放出了东方红的第一段，但是可以听出效果很不 好，只能听出具有东方红的调子而已。1.2 除噪音，加包络在1.1中听到有“啪”的杂声，下面通过加包络来消噪音。最简单的包络为指数衰减。最简单的指数衰减是对每个音乘以因子，在实验中首先加的是e5t的衰减，这种衰减方法使用的是相同速度的衰减，但是发现噪音并没有完全消 除，播放的音乐效果不是很好，感觉音乐起伏性不强。于是采用不同速度的衰减，根据乐音 持续时间的长短来确定衰减的快慢，乐音持

6、续时间越长，衰减的越慢，持续时间越短，衰减 的越快。在1.1程序的基础上加上包络，编写如下程序：soun d_1_21.mclear;clc;fs=8000;%由样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率time=fs*1/2,1/4,1/4,1,1/2,1/4,1/4,1;%各个乐音的抽样点数N=le ngth(time);%这段音乐的总抽样点数east=zeros(1,N);%用east向量来储存抽样点n=1;for num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/fs:1/fs:tim

7、e(num)/fs;%产生第nun个乐音的抽样点G=zeros(1,time( num);%为存储包络数据的向量G(1:time( nu m)=exp(1:(-1/time( nu m):1/8000);%产生包络点east( n:n+time( num)-1)=s in (2*pi*f( num)*t).*G(1:time( num);%合第nun个乐音加上包络n=n+time( nu m);endsou nd(east,8000);%播放plot(east);播放后可以听出噪音已经消除，同时因为不同时长的乐音衰减的快慢不一样，音乐听起 来更有起伏感，下图是加包络后的east图像。更科学的包

8、络如下图所示，每个乐音都经过冲激、衰减、持续、消失四个阶段。由上图可以看出这个包络是四段直线段构成的，因此只要确定了每段线段的端点，即可 用端点数据写出直线方程，因为直线方程可以用通式写出(我用的是斜截式)，因此这段包 络可以用简单的循环来完成。例如认为包络线上的数据如下图所示： 据此在MATLA中编写如下程序：sound_1_22.mclear;clc;fs=8000;%由样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;time=fs*1/2,1/4,1/4,1,1/2,1/4,1/4,1;N=le ngth(time);east

9、=zeros(1,N);n=1;for num=1:N%各个乐音对应的频率 %各个乐音的抽样点数%这段音乐的总抽样点数%用east向量来储存抽样点 %利用循环产生抽样数据，num表示乐音编号t=1/fs:1/fs:(time(num)/fs;%产生第 num个乐音的抽样点P=zeros(1,time( num);%为存储包络数据的向量L=(time( nu m)*0 1/5 333/1000 333/500 1;T=0 1.5 1 1 0;s=1;b=1:1:time( nu m);for k=1:4P(s:L(k+1)-1)=(T(k+1)-T(k)/(L(k+1)-L(k)*(b(s:L(

10、k+1)-1)-L(k+1)*o nes(1,L(k+1)-s)+T(k+1)*o nes(1,L(k+1)-s);%包络线直线方程通式沪生包络线抽样点%包络线端点对应的横坐标 %包络线端点对应的纵坐标s=L(k+1);end east( n:n+time( num)-1)=s in (2*pi*f( nu m)*t).*P(1:time( num);%合第nun个乐音加上包络n=n+time( nu m); endsou nd(east,8000); plot(east);运行得到的图像为：下图是两个乐音交接处的局部放大图，可以看到前一个乐音一直衰减到0,后一个乐音从0开始增加，因此消除了噪

11、音。若不需要每个音都衰减到0,例如只需衰减到持续阶段幅值的20%对程序做简单的修改 即可，将T=0 1.5 1 1 0改为T=0.2 1.5 1 1 0.2运行得到的结果为：由图可见，每个乐音都是衰减到一较小值而不是 0,也能消除噪音，同时音乐听起来更 加连续。1.3 改变程序，实现1.2中的音乐升高和降低一个八度升高一个八度即每个乐音的频率都提高一倍，变为原来的 2被；降低一个八度即每个乐 音的频率都减小一倍，变为原来的1/2。因此最简单的办法是将存储乐音频率的向量每个元 素改变为2或1/2倍。即将程序中的 f=523.25 523.25 587.33 392 349.23 349.23 2

12、93.66 392;改为f=523.25 523.25 587.33 392 349.23 349.23 293.66 392*2;或f=523.25 523.25 587.33 392 349.23 349.23 293.66 392/2;1/12将上述音乐上高半个音阶，即将频率变为原来的 2 (1.06 )倍，可以利用 resamlpe函数对原来的数据点进行重采样来实现east=resample(east,100,106);因为resample进行重新采样后会使每个乐音的持续时间改变，但是因为升高半个音阶，频率改变不大，所以每个音的持续时间是基本不变的。1.4 在 1.2 的音乐中加入谐波

13、在1.2 的音乐中加上二、三、四次谐波，基波幅度为 1，高次谐波幅度分别为 0.2 、0.3 、0.1 。只需将 1.2程序改为sound_1_4.mclear;clc;fs=8000;%由样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率 %各个乐音的抽样点数time=fs*1/2,1/4,1/4,1,1/2,1/4,1/4,1;N=le ngth(time);east=zeros(1,N);n=1;for num=1:N%这段音乐的总抽样点数%用east向量来储存抽样点%利用循环产生抽样数据，num表示乐音编号

14、t=1/fs:1/fs:(time(num)/fs;%产生第 num个乐音的抽样点P=zeros(1,time( num);%为存储包络数据的向量L=(time( nu m)*0 1/5 333/1000 333/500 1;%包络线端点对应的横坐标%包络线端点对应的纵坐标T=0 1.5 1 1 0;s=1;b=1:1:time( nu m);for k=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)*o nes(1,L(k+1)-s)+T(k+1)*o nes(1,L(k+1)-s);%包络线直线方程通式沪生包

15、络线抽样点s=L(k+1);endm=1 0.3 0.2; ss=zeros(1,le ngth(t);for i=1:length(m)ss=ss+m(i)*si n(2*i*pi*f( nu m)*t);endeast( n:n+time( nu m)-1)=ss.*P(1:time( nu m);%合第num个乐音加上包络澈形幅值矩阵%加谐波n=n+time( nu m); endsoun d(2*east,8000); plot(east);即可，加颜色部分为修改的部分，加上谐波后音乐效果变得更好了。1.5自选音乐合成两只老虎两只老虎曲调为C,因此可以得到每个乐音对应的频率分别为:12

16、345262.63293.66329.63349.23392每小节有四拍，一拍的时间是0.5s，因此各乐音的持续时间为:乐音12311231时间0.50.50.50.50.50.50.50.5乐音345345时间0.50.510.50.510.25sound 15.mclear;clc;fs=8000;%由样频率f=262.63 293.66 329.63 262.63 262.63 293.66 329.63 262.63 329.63349.23 392;%各个乐音对应的频率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,

17、1;%各个乐音的抽样点数N=le ngth(time);%这段音乐的总抽样点数east=zeros(1,N);%用east向量来储存抽样点n=1;for num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/fs:1/fs:time(num)/fs;%产生第nun个乐音的抽样点G=zeros(1,time( num);%为存储包络数据的向量G(1:time( nu m)=exp(1:(-1/time( nu m):1/8000);%产生包络点east( n:n+time( num)-1)=s in (2*pi*f( num)*t).*G(1:time( num);%合第nun个乐音加

18、上包络n=n+time( nu m);endsou nd(east,8000);%播放plot(east);第二部分用傅里叶变换分析音乐2.1 载入fmt.wav并播放利用wavread函数载入，用sou nd函数播放，程序如下：soun d 2 1.mwave=wavread('fmt.wav');soun d(wave)这段音乐听起来比之前合成的音乐更加真实，因为里边含有丰富的谐波。2.2 载入文件Guitar.mat，处理原始数据realwave载入文件Guitar.mat，分析wave2proc是怎么由realwave得到的。利用load Guitar.mat;载入并用

19、plot函数将realwave、wave2proc分别画出，得到以下两幅 图可以看到，wave2proc比realwave的周期性好得多，去掉了非线性谐波和噪声。在时域 做，从图上可以看到，realwave的数据大约是10个周期的共243个数据，因此可以用resample 函数对realwave进行重新采样，将采样点提高到250个，那么重采样后每个周期有25个点， 将这25个点对应相加求平均值后得到一个周期的值，因为进行了平均，减小了非线性谐波和 噪音，然后将这25个数据延托成十个周期即250个点，在利用resample函数对得到的函数重 新采样将采样点数恢复到243个。根据以上分析，编写实现

20、这个思路的程序如下：sound_2_2.mclear;clc;%重采样，将点数变为250%1个周期的对应点分别求和load Guitar.mat ;wave=resample(realwave,250,243); w=zeros(1,25);for i=1:25for k=0:9w(i)=w(i)+wave(25*k+i);endw=w/ 运仃后Iwave.wa%取平均值end,mat(w,1,10);%各1个周期的10个点延拓至250个点由图可见，两组组数据重合4的2很好;，说%重这样方将是很不错的方法。hold on,plot(wave2, 'r' ),hold off ;

21、 %各处理后的数据绘出，红色 h?ld on,plo分析wavc)2pro%将所的基波和谐波蓝色为了分析wave2proc的基波和谐波，可以对wave2proc进行傅里叶变换，得到wave2proc 的幅值谱，在频谱图上的第一个突出的波峰对应的频率即为wave2proc基频，利用help fft学习了 MATLA中快速傅里叶变换函数fft的用法，编写了如下程序：clear;clc;load Guitar.mat;fs=8000;NFFT = 25 extpow2(le ngth(wave2proc);Y = fft(wave2proc,NFFT)/le ngth(wave2proc);g =

22、fs/2*li nspace(0,1,NFFT/2+1); plot(g,2*abs(Y(1:NFFT/2+1)运行后得到的结果为虽然从图上可以大概看出包络，但是非常不明显，假如提高频域的抽样频率，例如将抽样频率由 NFFT = 2Anextpow2(length(wave2proc)改为NFFT = 8八门extpow2(length(wave2proc)得到的结果如下；由图可见虽然频域的抽样频率提高了很多，但是得到的包络依然不精确，这是因为 wave2proc是周期函数，但是现在的wave2proc只有243个数据点，并不能非常明显的体现出其周期性，因此它的幅值谱的离散化程度不高，虽然提高

23、了频域的抽样频率，但是wave2proc 数据点的周期性并没有增加，所以要显示出离散化程度高的幅值谱，就要增加wave2proc的周期性，即让wave2proc在时域重复多次后在进行傅里叶变换。利用repmat函数可以将wave2proc在时域重复。将程序修改为sound_2_3.mclear;clc;load Guitar.mat;fs=8000;wave2proc =repmat(wave2proc,20,1);%各 wave2proc重复 20次NFFT = 25 extpow2(le ngth(wave2proc);Y = fft(wave2proc,NFFT)/le ngth(wav

24、e2proc);g = fs/2*li nspace(0,1,NFFT/2+1);plot(g,2*abs(Y(1:NFFT/2+1)运行后得到的幅值谱为可以看出幅值谱的离散化程度已经非常高了。由图读出wave2proc的基频为329.1Hz，幅值为0.05401，高次谐波幅值分别为：谐 波23456789幅 值0.076760.048410.051900.0057090.019230.0067410.0073262.4 自动分析fmt.wav的音调和节拍思路分析：将fmt.wav导入后得到的是一个向量，它包含了这段音乐的所有信息，要自动分析这段 音乐的音调就需要将每个音调对应的点进行傅里叶变

25、换得到其幅值谱，在幅值谱上找到第一个幅值较大的极大值点，该点对应的就是该音调的基频，得到基频后就可以得到高次谐波的 幅值。为了使对每个音调进行傅里叶变换后得到的幅值谱离散程度高，应该将每个音调的数据在时域上重复多次，由于这些点都是直接采集的为做处理的点，因此其重复次数应该足够大才能体现出较强周期性，本实验采用重复1000次,虽然重复次数越多越好，但是次数太大， 程序运行的速度会大大降低。这里边还有两个关键点：第一，在从幅值谱上找基频时，因为图上的极大值点很多，怎 么能让程序自动确定出准确的基频。第二，在程序找到了基频之后，再由基频去获取高次谐 波的幅值时需要有一定的容错能力，例如若基频为 20

26、0Hz,幅值为1,那么对应的二次谐波的 频率为400Hz,但是很可能恰好幅值谱上400Hz处的幅值为0.01，但是401Hz处的幅值为0.2， 这时实际上的二次谐波应该为401Hz,但是若没有给基频一个容错范围，显然找到的二次谐 波的幅值是不正确的。针对以上提出的两个关键点，我找到了两条有针对性的解决办法。对于第一点，因为幅 值谱上极大值点的幅值足够大才能将其定位基频，因此在分析了几个音调后发现基频处的幅值都在0.025以上，因此将基频处的限定条件改为幅值大于 0.025的，但是在运行后发现，有 几个音调没有分析出来，说明它们的基频幅值小于0.025，其实可以观察一下fmt.wav的波形 就会

27、发现，有几段的整体幅值很小，因此基频幅值小，于是又在加上限定条件，若所有点的 幅值都小于0.02，那么再用0.015作为幅值的限定条件继续找，这样就将剩下的音调基频也 确定出来了。对于上述的第二点，可以将确定出的基频的误差设为+-1Hz，例如程序确定的基频为200Hz,实际的基频应该在(200-1 )到(200+1)之间，那么k次谐波对应的频率范围 是k*(200-1)到k*(200+1)，在这个区间中继续找幅值的极大值点就是k次谐波对应点。根据以上思路，下面开始编写用于分析一个音调频率的函数analysis 。在取谐波幅值时，幅值小于基波幅值5%勺谐波认为其幅值为0,最终谐波的幅值用归一 化

28、后的数据表示。每一步的详细思路见注释。an alysis.mfunction y1 y2=a nalysis(w,a)%5有两个返回值，y1返回频率，y2返回幅值，两个变量，w为待分析数组, a为数据的抽样频率fs=a;%傅里叶变换的抽样频率y仁zeros(1,7);%求最大7次谐波，因此定义1*7矩阵y2=zeros(1,7);NFFT = 2An extpow2(le ngth(w);Y = fft(w,NFFT)/le ngth(w);写完分manag = fs/2*li nspace(0,1,NFFT/2+1);Pld=首先time=W(Kn096；p)267 1.767 2.234

29、2.706 3.146 3.606.4.056 4P20>0I®2)&&(49P7810&(7(7)>P7.923L.%寻找基频的条件8.0281；1490(8)959 9.454 9.852 10%2储 1基频56.10.5y2(1)=p(k);11.292 11.741 12%存4储基波幅值 13.2693%； /16.3834*N '聊7若未找到基频，将幅'且限制改为 0.间转换成对应的数据点数，由于点的个数必须是整数，因此用 了对应音调)交接处的点后就可以进行编程了，用循环一个一个音调分析。程序如下sou nd_2_4.me

30、ndCeat=27wave=wa=teadd挪(綢);%卖入在N=le ngthfave);%确定数据总数(Pme:>1O(02)®&(P6i)>0651P767&&23(b>P(06)3&&(P&>O6+1)4.05筒嘶802费(4=0签9510.565 10.8g213.269商8 14.315 14.939 15.432/16.384*N); (p=N>10084&&(p(i)>(严買yg(脾(则2(1); for k=1yn(t)=g(i);if k=1高次谐波幅值temp=w

31、pve(1：time(k)-1);14.315 14.9391603S45是&&续寻找一项是为了将时 floor函数进行取整。在得到差允许的范围内寻找t次谐波点®454 9.85；幕站%!綁值归一化 eak.292 11.741 12.284 12.741.%节点向量0%5定数据的抽样频率 1)&&( p(i)>p(i+1)%这一段循环是与找基频条件相对应的获取运行后的结果见Excel文件“自动分析结果”。第三部分基于傅里叶级数的音乐合成3.1用2.3分析出来的结果重新加谐波基频329.1Hz 幅值为0.05401谐 波23456789幅 值0.

32、076760.048410.051900.0057090.019230.0067410.007326再次完成1.4只需将1.4程序中的波形幅度矩阵m=1 0.3 0.2改为m=0.05401 0.07676 0.04841 0.0519 0.005709 0.01923 0.006791 0.007326; 即可3.2 通过2.4提取的吉他音调信息弹奏东方红根据2.4分析的结果可以提取出吉他的音调信息，在“自动分析结果”Excel文件中用黄色标出了出来，一下所用的音调信息用最接近的频率来近似。而所给信息中缺少523.25和587.33的数据，因此用它们的一般来近似频率293.66349.23392523.25587.33近似频率291.97350.65391.94:2*261.632*295.4511111120.520305443004000050.05280271060.09594026000070000将1.4中的程序改为sound_3_2.mfs=8000;%由样频率f=523.25 523.25 587.33 392 349.23 349.23