Matlab音乐合成实验报告.docx

1、Matlab音乐合成实验报告音乐合成实验音乐合成实验 1摘要： 2第一部分 简单的合成音乐 21.1合成东方红 21.2 除噪音，加包络 41.3改变程序，实现1.2中的音乐升高和降低一个八度 111.4在1.2的音乐中加入谐波 121.5自选音乐合成一一两只老虎 1.3第二部分 用傅里叶变换分析音乐 1.5.2.1载入 fmt.wav 并播放 .15.2.2载入文件Guitar.mat，处理原始数据realwave 162.3分析wave2proc 的基波和谐波 1.92.4自动分析fmt.wav的音调和节拍 23第三部分 基于傅里叶级数的音乐合成 273.1用2.3分析出来的结果重新加谐波

2、 273.2 通过2.4提取的吉他音调信息弹奏东方红 2 8实验收获 31.摘要：本文共有三大部分：第一部分，简单的音乐合成；第二部分，用傅里叶变换 分析音乐；第三部分，基于傅里叶级数的音乐合成。由潜入深，一步一步分析了 用MATLAB进行音乐合成的过程。通过本实验达到了加深对傅里叶级数和傅里 叶分析的理解，熟悉对 MATLAB基本使用的目标。第一部分 简单的合成音乐1.1 合成东方红根据东方红第一小节的简谱和十二平均律计算出该小节每个乐音的频率， 在MATLAB中生成幅度为1，抽样频率为8kHz的正弦信号表示这些乐音，用 sou nd播放合成的音乐】=F十r 1 2 11邑i2 |4帽花3城

3、曲一巾务曲就由图可知东方红的曲调定为 F,即仁F,对应的频率为349.23Hz ,据此可以计算出其他乐音的频率，例如 5对应的频率为f5 349.23 27/12 523.25，一次类推计算出第一小节各乐音对应的频率为：乐音55621162频率523.2523.2587.3392349.2349.2293.6392* *553336在确定了各乐音的频率之后需要确定每个乐音的持续时间。每小节有两拍,一拍的时间是0.5s，因此各乐音的持续时间为：乐音55621162时间0.50.250.2510.50.250.251而在MATLAB中表示乐音所用的抽样频率为fs=8000Hz，也就是所1s钟内有8

4、000个点，抽样点数的多少就可表示出每个乐音的持续时间的长短。用一个行 向量来存储这段音乐对应的抽样点，在用 sou nd函数播放即可。根据以上分析在MATLAB中编写如下程序：soun d_1_1.mclear;clc;fs=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); east=zeros(1,N); n=1;腕段音乐的总抽样点数 %用 east向量来储存抽样点fo

5、r num=1:N%利用循环产生抽样数据，nun表示乐音编号t=1/fs: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);%播放音乐在MATLAB中运行sound_1_1.m，播放出了东方红的第一段，但是可以 听出效果很不好，只能听出具有东方红的调子而已。1.2 除噪音，加包络在1.1中听到有“啪”的杂声，下面通过加包络来消噪音。最简单的包络为指数衰减。最简单的指数衰减是对每个音乘以

6、e t因子，在1 51实验中首先加的是e .的衰减，这种衰减方法使用的是相同速度的衰减，但是 发现噪音并没有完全消除，播放的音乐效果不是很好，感觉音乐起伏性不强。于 是采用不同速度的衰减，根据乐音持续时间的长短来确定衰减的快慢，乐音持续 时间越长，衰减的越慢，持续时间越短，衰减的越快。在 1.1程序的基础上加上包络，编写如下程序：sound 1 21.mclear;clc;fs=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;

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

8、un个乐音加上包络n=n+time( nu m);endsou nd(east,8000); %播放plot(east);播放后可以听出噪音已经消除，同时因为不同时长的乐音衰减的快慢不一样,音乐听起来更有起伏感，下图是加包络后的 east图像。更科学的包络如下图所示，每个乐音都经过冲激、衰减、持续、消失四个阶由上图可以看出这个包络是四段直线段构成的， 因此只要确定了每段线段的 端点，即可用端点数据写出直线方程，因为直线方程可以用通式写出（我用的是 斜截式），因此这段包络可以用简单的循环来完成。例如认为包络线上的数据如 下图所示:据此在MATLAB中编写如下程序:sound_1_22.mclea

9、r;clc;fs=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); %这段音乐的总抽样点数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); %为存储包络数据的向

10、量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:4%包络线端点对应的横坐标%包络线端点对应的纵坐标泸生包络线抽样点P(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(n

11、um);%合第nun个乐音加上包络n=n+time( nu m);end sou nd(east,8000); plot(east);运行得到的图像为:21.510.50-0.5-1-1.5-20 0.5 1 1.5 2 2.5 3 3.54X 10下图是两个乐音交接处的局部放大图，可以看到前一个乐音一直衰减到 后一个乐音从0开始增加，因此消除了噪音。* *若不需要每个音都衰减到0，例如只需衰减到持续阶段幅值的20%，对程序做简单的修改即可，将T=0 1.5 1 1 0改为T=0.2 1.5 1 1 0.2运行得到的结果为：由图可见，每个乐音都是衰减到一较小值而不是 0，也能消除噪音，同时音乐

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

14、%各个乐音对应的频率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个乐音的抽样点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(

15、nu m);for k=1:4%包络线端点对应的横坐标%包络线端点对应的纵坐标泸生包络线抽样点P(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)-1

16、)=ss.*P(1:time( nu m);n=n+time( nu m); endsoun d(2*east,8000); plot(east);即可，加颜色部分为修改的部分，加上谐波后音乐效果变得更好了1.5 自选音乐合成一一两只老虎两只老虎1-C- 123,12311345H 345-w4曲调为C,因此可以得到每个乐音对应的频率分别为:12345262.63293.66329.63349.23392一拍的时间是因此各乐每小节有四拍，0.5s，音的持续时间为：乐音 12311231时间 0.50.50.50.50.50.50.50.5乐音 345345时间 0.50.510.50.510.

17、25sound_1_5.mclear;clc;fs=8000; %t样频率f=262.63 293.66 329.63 262.63 262.63 293.66 329.63 262.63 329.63 349.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/

18、fs:time( num)/fs; %产生第nun个乐音的抽样点G=zeros(1,time( num); %为存储包络数据的向量G(1:time( nu m)=exp(1:(-1/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函数播

19、放，程序如下soun d_2_1.mwave=wavread(fmt.wav);soun d(wave)这段音乐听起来比之前合成的音乐更加真实，因为里边含有丰富的谐波2.2 载入文件Guitar.mat ，处理原始数据 realwave载入文件Guitar.mat，分析wave2proc是怎么由realwave得到的。利用load Guitar.mat;载入并用 plot 函数将 realwave、wave2proc 分别画出，得到以下两幅图0.25wave2proc0.2i0.10.05II可以看到，wave2proc比realwave的周期性好得多，去掉了非线性谐波和 噪声。在时域做，从图

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

21、数变为 250w=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, r ),hold off ; %将处理后的数据绘出，红色 hold on,plot(wave2proc); %将所给的数据绘出，蓝色运行后的结果为:0.25由图可见，两组数据重合的很好，说明这种方

22、法是很不错的方法。2.3 分析wave2proc 的基波和谐波为了分析wave2proc的基波和谐波，可以对wave2proc进行傅里叶变换，得到wave2proc的幅值谱，在频谱图上的第一个突出的波峰对应的频率即为wave2proc基频，利用help fft学习了 MATLAB中快速傅里叶变换函数fft的用法,编写了如下程序:clear;clc;load Guitar.mat;fs=8000;NFFT = 25 extpow2(le ngth(wave2proc);Y = fft(wave2proc,NFFT)/le ngth(wave2proc); g = fs/2*li nspace(0

23、,1,NFFT/2+1);plot(g,2*abs(Y(1:NFFT/2+1)运行后得到的结果为0.08-b IIb n0.07-0.06-0.05-0.04-J111-0.031r 1-0.0211 I1I10.01卜/I ! |1 1 _7 J1 h1 |/ n Jr B1-J -100 500 1000 1500 2000 2500 3000 3500 4000虽然从图上可以大概看出包络，但是非常不明显，假如提高频域的抽样频率,例如将抽样频率由 NFFT = 2Anextpow2(length(wave2proc) 改为NFFT = 8八门extpow2(length(wave2proc

24、) 得到的结果如下；由图可见虽然频域的抽样频率提高了很多，但是得到的包络依然不精确，这 是因为wave2proc是周期函数，但是现在的wave2proc只有243个数据点，并 不能非常明显的体现出其周期性，因此它的幅值谱的离散化程度不高，虽然提高 了频域的抽样频率，但是wave2proc数据点的周期性并没有增加，所以要显示 出离散化程度高的幅值谱，就要增加 wave2proc的周期性，即让wave2proc在 时域重复多次后在进行傅里叶变换。利用repmat函数可以将wave2proc在时域重复。将程序修改为sound 2 3.mclear;clc;load Guitar.mat;fs=800

25、0;wave2proc =repmat(wave2proc,20,1); %各 wave2proc重复20次NFFT = 25 extpow2(le ngth(wave2proc);Y = fft(wave2proc,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.1Hz，幅值为0.05401，高次谐波幅值分别为:谐波23456789幅0.07670.04840.05100.005700.01

26、920.006740.00732值61993162.4 自动分析fmt.wav 的音调和节拍思路分析：将fmt.wav导入后得到的是一个向量，它包含了这段音乐的所有信息，要自动分析这段音乐的音调就需要将每个音调对应的点进行傅里叶变换得到其幅值谱，在幅值谱上找到第一个幅值较大的极大值点， 该点对应的就是该音调的基频， 得到基频后就可以得到高次谐波的幅值。为了使对每个音调进行傅里叶变换后得 到的幅值谱离散程度高，应该将每个音调的数据在时域上重复多次， 由于这些点都是直接采集的为做处理的点，因此其重复次数应该足够大才能体现出较强周期 性，本实验采用重复1000次，虽然重复次数越多越好，但是次数太大，

27、程序运 行的速度会大大降低。这里边还有两个关键点：第一，在从幅值谱上找基频时，因为图上的极大值 点很多，怎么能让程序自动确定出准确的基频。第二，在程序找到了基频之后， 再由基频去获取高次谐波的幅值时需要有一定的容错能力，例如若基频为200Hz，幅值为1，那么对应的二次谐波的频率为400Hz，但是很可能恰好幅值谱上 400Hz处的幅值为0.01，但是401Hz处的幅值为0.2，这时实际上的二次谐波应 该为401Hz，但是若没有给基频一个容错范围，显然找到的二次谐波的幅值是不 正确的。针对以上提出的两个关键点，我找到了两条有针对性的解决办法。对于第一 点，因为幅值谱上极大值点的幅值足够大才能将其定

28、位基频， 因此在分析了几个音调后发现基频处的幅值都在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，最终谐波的 幅值用归一化后的数据表示。每一步的详细思路见注释。an alysis.m fun cti on y1 y2=a nalysis(w,a) 做有两个返回值，y1返回频率，y2返回幅值，两个变量，w为待分析数组, a为数据的抽样频率fs=a; %f里叶变换的抽样频率y仁zeros(1,7); %求