Maab音乐合成实验报告

1、音乐合成实验目录音乐合成实验1摘要:1第一部分简单的合成音乐11.1 合成东方红 21.2 除噪音，加包络31.3 改变程序，实现1.2中的音乐升高和降低一个八度 71.4 在1.2的音乐中加入谐波81.5 自选音乐合成一两只老虎 8第二部分用傅里叶变换分析音乐92.1 载入fmt.wav并播放92.2 载入文件Guitar.mat,处理原始数据realwave 92.3 分析wave2proc的基波和谐波 112.4 自动分析fmt.wav的音调和节拍 14第三部分基于傅里叶级数的音乐合成163.1 用2.3分析出来的结果重新加谐波163.2 通过2.4提取的吉他音调信息弹奏东方红 17实验

2、收获17摘要：本文共有三大部分：第一部分，简单的音乐合成；第二部分，用傅里叶变换分析 音乐；第三部分，基于傅里叶级数的音乐合成。由潜入深，一步一步分析了用MATLA 进行音乐合成的过程。通过本实验达到了加深对傅里叶级数和傅里叶分析的理解，熟 悉对MATLA基本使用的目标。clear;clc;fs=8000; %t 样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率time=fs*1/251/431/4,1,1/2,1/4,1/4,1;% 各个乐音的抽样点数第一部分1.1 合成东方红简单的合成音乐根据东方红第一小节的

3、简谱和十二平均律计算出该小节每个乐音的频率，在 MATLAB八生成幅度为1才由样频率为8kHz的正弦信号表示这些乐音，用sound播放合 成的音乐- 5 Si 2 - If R | 2 - I4rfr-3磅曲（航方饪水$曲演由图可知东方红的曲调定为F,即仁F,对应的频率为349.23Hz,据此可以 计算出其他乐音的频率，例如5对应的频率为7/12fs 349.23 2523.25, 一次类推计算出第一小节各乐音对应的频率为：乐音55621162频率523.25523.25587.33392349.23349.23293.66392在确定了各乐音的频率之后需要确定每个乐音的持续时间。每小节有两拍

4、，一拍 的时间是0.5s,因此各乐音的持续时间为：乐音55621162时间0.50.250.2510.50.250.251而在MATLA中表示乐音所用的抽样频率为fs=8000Hz,也就是所1s钟内有 8000个点，抽样点数的多少就可表示出每个乐音的持续时间的长短。用一个行向量 来存储这段音乐对应的抽样点，在用sou nd函数播放即可。根据以上分析在MATLA中编写如下程序：soun d_1_1 ,mN=le ngth(time);%这段音乐的总抽样点数east=zeros(1,N);%用east向量来储存抽样点n=1；for num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/f

5、s:1/fs:time( num)/fs;%产生第nun个乐音的抽样点east( n:n+time( nu m)-1)=s in (2*pi*f( nu m)*t);%抽样点对应的幅值 n=n+time( nu m);endsou nd(east,8000); %播放音乐在MATLA中运行sound_1_1.m,播放出了东方红的第一段，但是可以听出 效果很不好，只能听出具有东方红的调子而已。1.2除噪音，加包络在1.1中听到有“啪”的杂声，下面通过加包络来消噪音最简单的包络为指数衰减。最简单的指数衰减是对每个音乘以 et因子，在实验中首先加的是的衰减，这种衰减方法使用的是相同速度的衰减，但是发

6、现 噪音并没有完全消除，播放的音乐效果不是很好，感觉音乐起伏性不强。于是采用 不同速度的衰减，根据乐音持续时间的长短来确定衰减的快慢，乐音持续时间越 长，衰减的越慢，持续时间越短，衰减的越快。在1.1程序的基础上加上包络，编写如下程序：clear;clc;fs=8000; %t 样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率time=fs*1/251/431/4,1,1/2,1/4,1/4,1;% 各个乐音的抽样点数sound 1 21.mclear;clc;fs=8000; %t 样频率f=523.25 52

7、3.25 587.33 392 349.23 349.23 293.66 392;%各个乐音对应的频率time=fs*1/251/431/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个乐音的抽样点G=zeros(1 ,time( num);%为存储包络数据的向量G(1 :time( nu m)=exp(1 :(-1/time( nu

8、m):1/8000);%产生包络点east( n:n+time( nu m)-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图像。210-1-24X 10更科学的包络如下图所示，每个乐音都经过冲激、衰减、持续、消失四个阶阕6. 5科*变化由上图可以看出这个包络是四段直线段构成的，因此只要确定了每段线段的 端点，即可

9、用端点数据写出直线方程，因为直线方程可以用通式写出（我用的是 斜 截式），因此这段包络可以用简单的循环来完成。例如认为包络线上的数据如下图 所示：据此在MATLA中编写如下程序: sound 1 22.mN=le ngth(time);%这段音乐的总抽样点数east=zeros(1,N);%用east向量来储存抽样点n=1；for num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/fs:1/fs:(time(num)/fs;%产生第 nun 个乐音的抽样点P=zeros(1 ,time( num);%为存储包络数据的向量L=(time( nu m)*0 1/5 333/1000

10、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 )* on es(1,L(k+1)-s);%包络线直线方程通式s=L(k+1);end east(n:n+time(num)-1 )=sin(2*pi*f(num)*t),*P(1 :time(num);%合第nun个乐音加上包络n=n+time

11、( nu m);end sou nd(east,8000); plot(east);运行得到的图像为:00.511.522.533.5X 100,下图是两个乐音交接处的局部放大图，可以看到前一个乐音一直衰减到 后一个乐音从。开始增加，因此消除了噪音。2若不需要每个音都衰减到0,例如只需衰减到持续阶段幅值的20%对程序做 简单的修改即可，将T=0 1.5 1 1 0改为T=0.2 1.5 1 1 0.2 运行得到的结果 为：00.511.522.533.54x 1021.510.50 -0.5-1-1.5-2由图可见，每个乐音都是衰减到一较小值而不是。，也能消除噪音，同时音乐听起来更加连续。1.

12、3 改变程序，实现1.2中的音乐升高和降低一个八度升高一个八度即每个乐音的频率都提高一倍，变为原来的2被；降低一个八度即 每个乐音的频率都减小一倍，变为原来的1/2o因此最简单的办法是将存储乐音频率 的向量每个元素改变为2或1/2倍。即将程序中的 f=523.25 523.25 587.33 392 349.23 349.23 293.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;将上述音乐上高半个音阶，

13、即将频率变为原来的 2仙2( 1.06)倍，可以利用resamlpe函数对原来的数据点进行重采样来实现 east=resample(east,100,106);因为resample进行重新采样后会使每个乐音的持续时间改变，但是因为升高 半个音阶，频率改变不大，所以每个音的持续时间是基本不变的。1.4 在1.2的音乐中加入谐波在1.2的音乐中加上二、三、四次谐波，基波幅度为1,高次谐波幅度分别为0.2、0.3、0.1 o只需将12程序改为sound 1 4,mclear;clc;fs=8000;样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66

14、392;%各个乐音对应的频率各个乐音的抽样点数time=fs*1/2,1/4,1/4J1/2,1/4,1/4,1;憾段音乐的总抽样点数%用east向量来储存抽样点利用循环产生抽 样数据，num表示乐音编号t=1/fs:1/fs:(time(num)/fs; %产生第 num 个乐音的抽样点%为存储包络数据的向量N=le ngth(time);east=zeros(1 ,N);n=1；for num=1 :NP=zeros(1,time( nu m);L=(time( nu m)*0 1/5 333/1000 333/500 1;%包络线端点对应的横坐标T=0 1.5 1 1 0;s=1；b=1

15、: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)* on es(1 ,L(k+1)-s);%包络线直线方程通式%包络线端点对应的纵坐标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( num)*t);%加谐波endeast( n:n+time( nu m

16、)-1)=ss.*P(1 :time( nu m);%合第num个乐音加上包络 n=n+time( nu m); endsoun d(2*east,8000); plot(east);即可，加颜色部分为修改的部分，加上谐波后音乐效果变得更好了1.5自选音乐合成一一两只老虎两只老虎1 = C-11231| 12311345-1 345-|A曲调为C,因此可以得到每个乐音对应的频率分别为:12345262.63293.66329.63349.23392clear;clc;fs=8000; %t 样频率修262.63 293.66 329.63 262.63 262.63 293.66 329.63

17、 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,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/

18、time( num):1/8000);%产生包络点east( n:n+time( nu m)-1)=s in (2*pi*f( num)*t).*G(1 :time( num);%合第nun个乐音加上包络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 载

19、入文件Guitar.mat ,处理原始数据realwave载入文件Gu让ar.mat,分析wave2proc是怎么由realwave得到的。利用loadGuitar.mat;载入并用plot函数将realwave、wave2proc分别画出，得到以下两幅图0.250.20.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.20501001502002501 Hi-100150502002500.20.150.10.05-0.05-n 10.25-0.15 0可以看到，wave2proc比realwave的周期性好得多，去掉了非线性谐波和噪 声。 在时域做，从图上可以看到，re

20、alwave的数据大约是10个周期的共243个数据，因此 可以用resample函数对realwave进行重新采样，将采样点提高到250个，那么重采样 后每个周期有25个点，将这25个点对应相加求平均值后得到一个周期 的值，因为进 行了平均，减小了非线性谐波和噪音，然后将这25个数据延托成十个周期即250个点， 在利用resample函数对得到的函数重新采样将采样点数恢复到243个。根据以上分析, 编写实现这个思路的程序如下：sound_2_2.mclear;clc;load Guitar.mat ; wave=resample(realwave,250,243); %重采样，将点数变为 25

21、0 w=zeros(1,25);for i=1:25for k=0:9w(i)=w(i)+wave(25*k+i);%1个周期的对应点分别求和 endendw=w/10;%取平均值wave2=repmat(w, 1,10);%将1个周期的10个点延拓至250个点wave2=resample(wave2,243,250); %重采样，将点数变回 243hold on,plot(wave2, Y* ),hold off ; % 将处理后的数据绘出，红色 holdon,plot(wave2proc);%将所给的数据绘出，蓝色运行后的结果为:0.250.20.150.10.05 0-0.05-0.1-

22、0.15由图可见，两组数据重合的很好，说明这种方法是很不错的方法。2.3分析wave2proc的基波和谐波为了分析wave2proc的基波和谐波，可以对wave2proc进行傅里叶变换，得到 wave2proc的幅值谱，在频谱图上的第一个突出的波峰对应的频率即为wave2proc 基频，利用help fft学习了 MATLA中快速傅里叶变换函数fft的用法，编写了如 下程 序：clearjclc;load Guitar.mat;fs=8000;NFFT = 25 extpow2(le ngth(wave2proc);Y = fft(wave2proc,NFFT)/le ngth(wave2pr

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

24、离散化程度不高，虽然提高了频域 的抽样频率，但是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(wave2pro

25、c,NFFT)/le ngth(wave2proc);g = fs/2*li nspace(0,1 ,NFFT/2+1);plot(g,2*abs(Y(1 :NFFT/2+1)可以看出幅值谱的离散化程度已经非常高了。由图读出wave2proc的基频为329.1 H乙幅值为0.05401,高次谐波幅值分别为：3456幅 0.07676 0.04841值0.0073260.05190.005709 0.01923 0.00674102.4 自动分析fmt.wav的音调和节拍思路分析：将fmt.wav导入后得到的是一个向量，它包含了这段音乐的所有信息，要自动分 析这段音乐的音调就需要将每个音调对应的

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

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

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

29、 幅值用归一化后的数据表示。每一步的详细思路见注释。analysis.mfun cti on y1 y2=a nalysis(w,a)做有两个返回值，y1返回频率，y2返回幅值，两个变量，w为待分析数组，a为数据的抽样频率fs=a;%f里叶变换的抽样频率yfczeros(1,7); %求最大7次谐波，因此定义1/7矩阵y2=zeros(1,7);NFFT = 2An extpow2(le ngth(w);Y = fft(w,NFFT)/le ngth(w); g = fs/2*li nspace(0,1 ,NFFT/2+1);p=2*abs(Y(1:NFFT/2+1);plot(g,p)%以上

30、为傅里叶变换部分d=floor(NFFT/fs);%将误差1 Hz化成对应的点数for k=2:length(p)-1if (p(k)>0.02)&&(p(k)>p(k-1 )&&(p(k)>p(k+1)%寻找基频的条件y1(1)=g(k);烷储基频y26=p(k)；烷储基波幅值breakelseif (p(k)>0.015)&&(p(k)>p(k-1 )&&( p(k)>p(k+1)%若未找到基频，将幅值限制改为0.015,继续寻找yi(i)=g(k)；y2(1)=p(k);endendfor

31、 t=2:7for i=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); %谐波 幅值归一化yi(t)=g(i);breakelseif(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)

32、;y1 (t)=g(i);%这一段循环是与找基频条件相对应的获取高次谐波幅值endendendfun cti on y1 y2=a nalysis(w,a)做有两个返回值，y1返回频率，y2返回幅值，两个变量，w为待分析数组，a为数据的抽样频率fs=a;%f里叶变换的抽样频率y仁zeros。,7);%求最大7次谐波，因此定义1*7矩阵在编写完分析函数analysis后即可编写自动分析的主程序。首先在Coool Edit中手动标定音调交界处的时间节点，得到了 time向量time=floor(0.096 0.267 1.767 2.234 2.706 3.146 3.606.4.056 4.52

33、0 5.030 5.749 5.978 7.015 7.709 7.923.8.028 8.490 8.959 9.454 9.852 10.125 10.356.10.565 10.822 11.292 11.741 12.284 12.741.13.269 13.758 14.315 14.939 15.432/16.384*N);716.384*N '中16.384是fmt.wav的总长度，N为fmt.wav数据点的总数，这一项 是为了将时间转换成对应的数据点数，由于点的个数必须是整数，因此用 floor 函数进行取整。在得到了对应音调交接处的点后就可以进行编程了，用循环一个一个

34、音调分析。程序如下： sound_2_4.mclear;clc;wave=wavread (fmt.wav1);%读入文件N=le ngth(wave);%确定数据总数time=floor(0.096 0.267 1.767 2.234 2.706 3.146 3.606.4.056 4.520 5.030 5.749 5.978 7.015 7.709 7.923.8.028 8.490 8.959 9.454 9.852 10.125 10.356.10.565 10.822 11.292 11.741 12.284 12.741.13.269 13.758 14.315 14.939 1

35、5.432/16.384*N);%节点向量fs=N/16.384;%确定数据的抽样频率n=len gth(time);for k=1 :n; if k=1temp=wave(1 :time(k)-1);else temp=wave(time(k-1):time(k)-1);%各第k个音调数据存入temp矩阵 end temp=repmat(temp,1000,1);%各数据重复 1000 次F(k,1:7) U(k,1:7)=a nalysis(temp,fs);%各每个音调的处理结果分别保存，F保存频率，U保存幅值 end运行后的结果见Excel文件“自动分析结果”。第三部分 基于傅里叶级数

36、的音乐合成3.1 用2.3分析出来的结果重新加谐波基频329.1Hz幅值为0.05401幅 0.07676 0.04841 0.05190.005709 0.01923 0.006741 0.0073260值再次完成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文件中 用黄色标出了出来，一下所用的音调信息用最接近的频率

37、来近似。而所给 信息中缺 少523.25和587.33的数据，因此用它们的一般来近似频率293.66349.23392523.25587.33近似频率291.97350.65391.942*261.632*295.4511111120.520305440.2757865390.1086165950.3071496930.205838324300.1758122440.06907235600.15708081340000.057298602050.0528027110.05603520200.0753863870.08891702260.095940260000700000.065977510将1.4中的程序改为sound 3 2.mfs=8000; %t 样频率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);%这段音乐的总抽样