Matlab音乐合成试验报告

1、音乐合成实验目录音乐合成实验. 1摘要：. 1第一部分 简单的合成音乐. 21.1合成东方红. 21.2除噪音，加包络. 31.3改变程序，实现 1.2中的音乐升高和降低一个八度 . 81.4在 1.2的音乐中加入谐波. 81.5自选音乐合成- 两只老虎 . 9第二部分 用傅里叶变换分析音乐. 102.1载入 fmt.wav并播放 . 102.2载入文件 Guitar.mat，处理原始数据 realwave . 112.3分析 wave2proc的基波和谐波. 142.4自动分析 fmt.wav的音调和节拍 . 17第三部分基丁傅里叶级数的音乐合成. 203.1用 2.3分析出来的结果重新加谐

2、波 . 203.2通过 2.4提取的吉他音调信息弹奏东方红 . 21实验收获. 22摘要:本文共有三大部分：第一部分，简单的音乐合成；第二部分，用傅里叶变换 分析音乐； 第三部分， 基丁傅里叶级数的音乐合成。 由潜入深， 一步一步分析了 用 MATLAffi行音乐合成的过程。通过本实验达到了加深对傅里叶级数和傅里叶 分析的理解，熟悉对 MATLABS本使用的目标。第一部分简单的合成音乐1.1合成东方红根据东方红第一小节的简谱和十二平均律计算出该小节每个乐音的频率， 在MATLAB 生成幅度为 1,抽样频率为 8kHz 的正弦信号表示这些乐音，用 sound 播放合成的音乐5| 2 1 161

3、2 - |压曲（嘉床如:第小市曲谓由图可知东方红的曲调定为 F,即 1=F,对应的频率为 349.23Hz,据此可以计算出其他乐音的频率，例如 5对应的频率为7/12f5= 349.23 x 2= 523.25 ,一次类推计算出第一小节各乐音对应的频率为：乐音55621162频率523.25523.25587.33392349.23349.23293.66392在确定了各乐首的频率之后需要确定每个乐首的持续时间。每小节有两拍,一拍的时间是 0.5s ,因此各乐音的持续时间为:乐音55621162时间0.50.250.2510.50.250.251而在 MATLAB 表示乐音所用的抽样频率为 f

4、s=8000Hz,也就是所 1s 钟内有 8000个点，抽样点数的多少就可表示出每个乐音的持续时间的长短。用一个行向 量来存储这段音乐对应的抽样点，在用 sound 函数播放即可。根据以上分析在 MATLAB 编写如下程序：sound_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=length(time);%段音乐的总抽样点数east=zeros(1,N);咖 ea

5、st向量来储存抽样点n=1;for num=1:N呀 U 用循环产生抽样数据，nunmt 示乐音编号t=1/fs:1/fs:time(num)/fs;跑生第 num乐音的抽样点east(n:n+time(num)-1)=sin(2*pi*f(num)*t);盟 t样点对应的幅值n=n+time(num);endsound(east,8000);%8 放音乐在 MATLAB 运行 sound_1_1.m,播放出了东方红的第一段，但是可以听 出效果很不好，只能听出具有东方红的调子而已。1.2除噪音，加包络在 1.1中听到有“啪”的杂声，下面通过加包络来消噪音最简单的包络为指数衰减。最简单的指数衰减

6、是对每个音乘以e因子，在实验中首先加的是e.5t的衰减，这种衰减方法使用的是相同速度的衰减，但是发现噪音并没有完全消除，播放的音乐效果不是很好，感觉音乐起伏性不强。丁 是采用不同速度的衰减，根据乐音持续时间的长短来确定衰减的快慢，乐音持续 时间越长，衰减的越慢，持续时间越短，衰减的越快。在 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

7、,1;以各个乐音的抽样点数N=length(time);%段音乐的总抽样点数east=zeros(1,N);咖 east向量来储存抽样点n=1;for num=1:N呀 U用循环产生抽样数据，nunmt 示乐音编号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+ti

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

9、ound_1_22.mclear;clc;fs=8000;盟 t样频率f=523.25 523.25 587.33 392 349.23 349.23 293.66 392;潞个乐音对应的频率潞个乐音的抽样点数呀 U用循环产生抽样数据，numg 示乐音编号%产生第 num乐音的抽样点%的存储包络数据的向量T=0 1.5 1 1 0;s=1;b=1:1:time(num);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)*ones(1,L(k+1)-s)+T(k+1)*ones(1,L(k+1)-s);%

10、包络线直线方程通式s=L(k+1);end east(n:n+time(num)-1)=sin(2*pi*f(num)*t).*P(1:time(num);%合第 num乐音加上包络n=n+time(num); endsound(east,8000);plot(east);运行得到的图像为:time=fs*1/2,1/4,1/4,1,1/2,1/4,1/4,1;N=length(time);east=zeros(1,N);n=1;for num=1:N脂段音乐的总抽样点数 咖 east向量来储存抽样点t=1/fs:1/fs:(time(num)/fs;P=zeros(1,time(num);L

11、=(time(num)*0 1/5 333/1000 333/500 1;%包络线端点对应的横坐标%包络线端点对应的纵坐标%产生包络线抽样点4X 10下图是两个乐音交接处的局部放大图，可以看到前一个乐音一直衰减到0,后一个乐音从 0开始增加，因此消除了噪音。1.510.50-0.5-1-1.5-2- *1.522.53.5若不需要每个音都衰减到0,例如只需衰减到持续阶段幅值的20%对程序做 简单的修改即可，将 T=0 1.5 1 1 0 改为 T=0.2 1.5 1 1 0.2 运行得到的结果 为：21.510.50-0.5-1-1.5-2由图可见，每个乐音都是衰减到一较小值而不是0,也能消除

12、噪音，同时音乐听起来更加连续。00.511.522.533.54X 101.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.25523.25 587.33 392 349.23 349.23 293.66 392*2; 或f=523.25 523.25 587.33

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

14、 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=length(time);%段音乐的总抽样点数east=zeros(1,N);咖 east 向量来储存抽样点n=1;for num=1:N%用循环产生抽样数据，nunmt 示乐音编号t=1/fs:1/fs:(time(num)/fs;以产生第 num乐音的抽样点P=zeros(1,time(num);%的存储包络数据的向量L=(time(num)*0 1/5 333/1000 333/500

15、1;%包络线端点对应的横坐标T=0 1.5 1 1 0;%包络线端点对应的纵坐标s=1;b=1:1:time(num);%产生包络线抽样点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)*ones(1,L(k+1)-s)+T(k+1)*ones(1,L(k+1)-s);%包络线直线方程通式s=L(k+1); endm=1 0.3 0.2;瞰形幅值矩阵ss=zeros(1,length(t);for i=1:length(m)ss=ss+m(i)*sin(2*i*pi*f(num)*t);%J 口谐波en

16、d east(n:n+time(num)-1)=ss.*P(1:time(num); %合第 num乐音加上包络n=n+time(num);end sound(2*east,8000); plot(east);即可，加颜色部分为修改的部分，加上谐波后音乐效果变得更好了1.4在1.2的音乐中加入谐波两只老虎41 = C|1231|L231 34A|345-|，曲调为 C,因此可以得到每个乐音对应的频率分别为:12345262.63293.66329.63349.23392每小节有四拍，一拍的时间是 0.5s,因此各乐音的持续时间为:乐音12311231时间0.50.50.50.50.50.50.

17、50.5乐音345345时间0.50.510.50.510.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.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=length(time);%段音乐的总抽样点数east=zeros(1,N);咖 east 向量来储存抽样点n=1;for num=1:N呀 U用循环

18、产生抽样数据，nunmt 示乐音编号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);endsound(east,8000);%8放1.5自选音乐合成一一两只老虎plot(east);第二部分用傅里叶变换分析音乐2.1载入fmt.wav并播放利用

19、 wavread函数载入，用 sound函数播放，程序如下:sound_2_1.mwave=wavread(fmt.wav);sound(wave)这段音乐听起来比之前合成的音乐更加真实，因为里边含有丰富的谐波。2.2载入文件Guitar.mat，处理原始数据realwave载入文件 Guitar.mat，分析 wave2proc是怎么由 realwave 得到的。利用load Guitar.mat; 载入并用 plot 函数将 realwave、wave2proc分别画出，得 到以下两幅图0.250.20.150.10.050-0.05-0.1-0.15-0.2 0realwave50100

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

21、，编写实现这个思路的程序如下：sound_2_2.mclear;clc; load Guitar.mat ; wave=resample(realwave,250,243);咖采样，将点数变为250w=zeros(1,25); for i=1:25 for k=0:9 w(i)=w(i)+wave(25*k+i); %101周期的对应点分别求和end end w=w/10;%R 平均值wave2=repmat(w,1,10);%务 1 个周期的 10个点延拓至 250 个点wave2=resample(wave2,243,250); 咖采样，将点数变回 243 hold on,plot(wav

22、e2, r ),holdoff ; %务处理后的数据绘出，红色 hold on,plot(wave2proc);%务所给的数据绘出，蓝色0.20.150.10.050-0.05-0.150100150200250运行后的结果为:2.3分析wave2proc的基波和谐波为了分析 wave2proc的基波和谐波，可以对 wave2proc进行傅里叶变换，得到wave2proc的幅值谱，在频谱图上的第一个突出的波峰对应的频率即为wave2proc基频，利用 help fft 学习了 MATLAB 快速傅里叶变换函数 fft的用法，编写了如clear;clc;load Guitar.mat;0.250

23、.20.150.10.050-0.05-0.1-0.15由图可见，两组数据重合的很好，说明这种方法是很不错的方法。运行后得到的结果为fs=8000;NFFT = 2Anextpow2(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)运行后得到的结果为0.080.070.060.050.040.030.020.010运行后得到的结果为虽然从图上可以大概看出包络，但是非常不明显，假如提高频域的抽样频率,05001

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

25、lc;load Guitar.mat;fs=8000;wave2proc =repmat(wave2proc,20,1);%务 wave2proc重复 20 次NFFT = 2Anextpow2(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,高次谐波幅值分别为：谐波23456789幅值0.076

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

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

28、法。对丁第一 点，因为幅值谱上极大值点的幅值足够大才能将其定位基频，因此在分析了几个音调后发现基频处的幅值都在 0.025以上，因此将基频处的限定条件改为幅值大 丁0.025的，但是在运行后发现，有几个音调没有分析出来，说明它们的基频幅 值小丁0.025,其实可以观察一下 fmt.wav 的波形就会发现，有几段的整体幅值很 小，因此基频幅值小，丁是乂在加上限定条件，若所有点的幅值都小丁0.02,那么再用 0.015作为幅值的限定条件继续找，这样就将剩下的音调基频也确定出来 了。对丁上述的第二点，可以将确定出的基频的误差设为+-1Hz,例如程序确定的基频为 200H乙实际的基频应该在(200-1

29、 )至 U(200+1)之间，那么 k次谐波对 应的频率范围是 k*(200-1)到 k*(200+1)，在这个区间中继续找幅值的极大值点就 是 k次谐波对应点。根据以上思路，下面开始编写用丁分析一个音调频率的函数analysis。在取谐波幅值时，幅值小丁基波幅值 5%勺谐波认为其幅值为 0,最终谐波的 幅值用归一化后的数据表示。每一步的详细思路见注释。analysis.mfunction y1 y2=analysis(w,a)以投有两个返回值，y1返回频率,y2返回幅值，两个变量,w为待分析数组, a 为数据的抽样频率fs=a;%f里叶变换的抽样频率y1=zeros(1,7);咻最大 7 次

30、谐波，因此定义 1*7矩阵y2=zeros(1,7); NFFT = 2Anextpow2(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);%务误差 1H也成对应的点数for k=2:length(p)-1if (p(k)0.02)&(p(k)p(k-1)&(p(k)p(k+1)%!找基频的条件y1(1)=g(k);帝储基频y2(1)=p(k);帝储基波幅值break

31、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 for t=2:7for i=t*(k-d):t*(k+d)9 疑误差允许的范围内寻找 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)

32、&(p(i)p(i-1)&(p(i)p(i+1) y2(t)=p(i)/y2(1);y1(t)=g(i);腺一段循环是与找基频条件相对应的获取高次谐波幅值 end endend 在编写完分析函数 analysis 后即可编写自动分析的主程序。 首先在 Coool Edit 中手动标定音调交界处的时间节点，得到了 time 向量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.0157.709 7.923. 8.028 8.490 8.959 9.454 9.

33、852 10.125 10.356. 10.565 10.822 11.292 11.74112.284 12.741.13.269 13.758 14.315 14.939 15.432/16.384*N);/16.384*N 中 16.384 是 fmt.wav 的总长度，N 为 fmt.wav 数据点的总数，这一项 是为了将时间转换成对应的数据点数，由丁点的个数必须是整数，因此用 floor 函数进行取整。在得到了对应音调交接处的点后就可以进行编程了，用循环一个一个音调分析。程序如下： sound_2_4.mclear;clc;wave=wavread(fmt.wav);咐入文件N=le

34、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 15.432/16.384*N);%节点向量fs=N/16.384;%!定数据的抽样频率n=length(time);fork=1:n

35、;ifk=1temp=wave(1:time(k)-1);elsetemp=wave(time(k-1):time(k)-1);%务第 k 个音调数据存入 temp 矩阵 endtemp=repmat(temp,1000,1);%务数据重复 1000次F(k,1:7) U(k,1:7)=analysis(temp,fs);%务每个音调的处理结果分别保存,F 保存频率,U保存幅值 end运行后的结果见Excel文件自动分析结果”第三部分 基于傅里叶级数的音乐合成3.1用2.3分析出来的结果重新加谐波基频 329.1Hz 幅值为 0.05401谐波23456789幅值0.076760.048410

36、.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 文件中用黄色标出了出来，一下所用的音调信息用最接近的频率来近似。 而所给 信息中缺少5 2 3 . 2 5和5 8 7 . 3 3的 数 据 ， 因 此 用 它 们 的 一 般 来

37、近 似频率293.66349.23392523.25587.33近似频 率291.97350.65391.942*261.632*295.451 111r 1 1120.520305440.2757865390.1086165950.3071496930.205838324300.1758122440.06907235600.15708081340000.057298602050.052802710.05603520200.0753863870.0889170226:0.095940260r 0 100700000.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=length(ti