完整版使用matlab进行简单音乐合成

1、信号与系统-综合实验之音乐合成(1)请根据?东方红?片断的简谱和“十二平均律计算出该片断中各个乐音的频率,在MATLAB中生成幅度为1、抽样频率为8kHz的正弦信号表示这些乐音.请用sound函数播放每个乐音,听一听音调是否正确.最后用这一系列乐音信号拼出?东方红?片断,注意限制每个乐音持续的时间要符合节拍,用sound播放你合成的音乐,听起来感觉如何代码如下：f =8000;t2=0:1/f:1;t4=0:1/f:0.5;t8=0:1/f:0.25;omg5=523.35;omg6=587.33;omg2=392;omg1=349.23;omg6l=293.66;m1=sin(2*pi*om

2、g5*t4);m2=sin(2*pi*omg5*t8);m3=sin(2*pi*omg6*t8);m4=sin(2*pi*omg2*t2);m6=sin(2*pi*omg1*t4);m7=sin(2*pi*omg1*t8);m8=sin(2*pi*omg6l*t8);m9=sin(2*pi*omg2*t2);m=m1 m2 m3 m4 m6 m7 m8 m9;sound(m);听的时候发现在相邻乐音之间有杂音,这是由于相位不连续造成的.你一定注意到(1)的乐曲中相邻乐音之间有“啪的杂声,这是由于相位不连续 产生了高频分量.这种噪声严重影响合成音乐的质量,丧失真实感.为了消 除它,我们可以用图1

3、.5所示包络修正每个乐音,以保证在乐音的邻接处信号幅度为零.此外建议用指数衰减的包络来表示我采用的是指数衰减的包络. 代码如下：f =8000;t2=0:1/f:1;t4=0:1/f:0.5;t8=0:1/f:0.25;omg5=523.35;omg6=587.33;omg2=392;omg1=349.23;omg6l=293.66;m1=exp(-2*t4).*sin(2*pi*omg5*t4); m2=exp(-4*t8).*sin(2*pi*omg5*t8); m3=exp(-4*t8).*sin(2*pi*omg6*t8); m4=exp(-1*t2).*sin(2*pi*omg2*t

4、2); m6=exp(-2*t4).*sin(2*pi*omg1*t4); m7=exp(-4*t8).*sin(2*pi*omg1*t8); m8=exp(-4*t8).*sin(2*pi*omg6l*t8); m9=exp(-1*t2).*sin(2*pi*omg2*t2); m=m1 m2 m3 m4 m6 m7 m8 m9; sound(m);第一次我采用的指数衰减没有时间前面的系数,即每个都只乘expt,没有系数;后来根据不同节拍,更改了不同的衰减系数,这样声音听起来感觉更加圆润.3请用最简单的方法将2中的音乐分别升高和降低一个八度.提示：音乐播放的时间可以变化再难一些,请用 res

5、ample函数也可以用interp和decimate函数将上述音乐升高半个音阶.提示：视计算复杂度,不必特别精确答：最简单的方法是直接更改抽样频率f o将f从8徵为4K,那么升高一个八度,并 且播放速度增快了一倍；将f从8k改为16k,那么降低一个八度,速度也变慢了一倍.改为4000或者16000升高半个音阶,只须在2代码最后加一句resample m,1000,1059即可.f =8000;t2=0:1/f:1;t4=0:1/f:0.5;t8=0:1/f:0.25;omg5=523.35;omg6=587.33;omg2=392;omg1=349.23; omg6l=293.66;m1=ex

6、p(-2*t4).*sin(2*pi*omg5*t4);m2=exp(-4*t8).*sin(2*pi*omg5*t8);m3=exp(-4*t8).*sin(2*pi*omg6*t8);m4=exp(-1*t2).*sin(2*pi*omg2*t2);m6=exp(-2*t4).*sin(2*pi*omg1*t4);m7=exp(-4*t8).*sin(2*pi*omg1*t8);m8=exp(-4*t8).*sin(2*pi*omg6l*t8);m9=exp(-1*t2).*sin(2*pi*omg2*t2);m=m1 m2 m3 m4 m6 m7 m8 m9;resample(m,100

7、0,1059);sound(m);(4)试着在(2)的音乐中增加一些谐波分量,听一听音乐是否更有“厚度 了注意谐波分量的能量要小,否那么掩盖住基音反而听不清音调了.(如果选择基波幅度为1 ,二次谐波幅度0:2 ,三次谐波幅度0:3 ,听起来像不像象风琴) 代码如下： f =8000;t2=0:1/f:1;t4=0:1/f:0.5;t8=0:1/f:0.25;omg5=523.35;omg6=587.33;omg2=392;omg1=349.23;omg6l=293.66;m1=exp(-2*t4).*sin(2*pi*omg5*t4)+0.2*exp(-2*t4).*sin(2*pi*2*om

8、g5*t4)+0.3*exp(-2*t4).*sin(2*pi*3*omg5*t4);m2=exp(-4*t8).*sin(2*pi*omg5*t8)+0.2*exp(-4*t8).*sin(2*pi*2*omg5*t8)+ 0.3*exp(-4*t8).*sin(2*pi*3*omg5*t8);m3=exp(-4*t8).*sin(2*pi*omg6*t8)+0.2*exp(-4*t8).*sin(2*pi*2*omg6*t8)+ 0.3*exp(-4*t8).*sin(2*pi*3*omg6*t8);m4=exp(-1*t2).*sin(2*pi*omg2*t2)+0.2*exp(-1*t

9、2).*sin(2*pi*2*omg2*t2)+ 0.3*exp(-1*t2).*sin(2*pi*3*omg2*t2);m6=exp(-2*t4).*sin(2*pi*omg1*t4)+0.2*exp(-2*t4).*sin(2*pi*2*omg1*t4)+0.3*exp(-2*t4).*sin(2*pi*3*omg1*t4);m7=exp(-4*t8).*sin(2*pi*omg1*t8)+0.2*exp(-4*t8).*sin(2*pi*2*omg1*t8)+ 0.3*exp(-4*t8).*sin(2*pi*3*omg1*t8);m8=exp(-4*t8).*sin(2*pi*omg6

10、l*t8)+0.2*exp(-4*t8).*sin(2*pi*2*omg6l*t8 )+0.3*exp(-4*t8).*sin(2*pi*3*omg6l*t8);m9=exp(-1*t2).*sin(2*pi*omg2*t2)+0.2*exp(-1*t2).*sin(2*pi*2*omg2*t2)+0.3*exp(-1*t2).*sin(2*pi*3*omg2*t2);m=m1 m2 m3 m4 m6 m7 m8 m9;sound(m);参加谐波分量后,音色有所变化,感觉更加清脆一些.(5)自选其它音乐合成,例如贝多芬第五交响乐的开头两小节.我选取的是?晴天?的第一句代码如下：f=8000;t

11、2=0:1/f:1;t4=0:1/f:0.5;t8=0:1/f:0.25;t=0:1/f:0.125;omg1=392;omg2=440;omg3=493.88;omg4=523.25;omg5=587.33;omg6=659.25;omg7=698.45;omg5l=293.66;m0=0;m1=exp(-2*t4).*sin(2*pi*omg5*t4);m2=exp(-2*t4).*sin(2*pi*omg5*t4);m3=exp(-2*t4).*sin(2*pi*omg1*t4);m4=exp(-1*t2).*sin(2*pi*omg1*t2);m5=exp(-2*t4).*sin(2*

12、pi*omg2*t4);m6=exp(-2*t2).*sin(2*pi*omg3*t2);m7=exp(-2*t4).*sin(2*pi*omg5*t4);m8=exp(-2*t4).*sin(2*pi*omg5*t4);m9=exp(-2*t4).*sin(2*pi*omg1*t4);m10=exp(-2*t4).*sin(2*pi*omg1*t4);m11=exp(-4*t8).*sin(2*pi*omg2*t8);m12=exp(-4*t8).*sin(2*pi*omg3*t8);m13=exp(-4*t8).*sin(2*pi*omg2*t8);m14=exp(-4*t8).*sin(

13、2*pi*omg1*t8);m15=exp(-2*t2).*sin(2*pi*omg5l*t2);m=m0 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15;sound(m);(6)先用wavread函数载入光盘中的fmt.wav文件,播放出来听听效果如何是否比刚刚的合成音乐真实多了 x=wavread( 'fmt.wav' );sound(x);(7)你知道待处理的wave2proc是如何从真实值realwave中得到的么这个预处理过程可以去除真实乐曲中的非线性谐波和噪声,对于正确分析音调是非常重要的.提示：从时域做,可以继

14、续使用 resample函数.realwave中的波形有十个周期,要除去其中的噪声可以采用时域求均值的方 法.步骤如下：首先用resample函数将其采样率增大为十倍,再等分十份取平均, 然后重复这个平均后的波形十次,复原其长度,最后在用resample函数复原采样率到1/10.代码如下：load( 'Guitar.MAT' );wavetemp = zeros(length(realwave), 1);waveresampled = resample(realwave,10,1);for n = 1 : 10wavetemp = wavetemp + waveresample

15、d(n - 1) * length(realwave) + 1 : n* length(realwave) / 10;endmywave2proc = repmat(wavetemp, 10, 1);mywave2proc = resample(mywave2proc, 1, 10);figure;subplot(3, 1, 1);plot(realwave);subplot(3, 1,2);plot(wave2proc);subplot(3, 1, 3);plot(mywave2proc);0.5 11c1:Ar,iRftrikAft0 - aj 3小h V/j 小 -0.5 1:05010

16、0150200250(8)这段音乐的基频是多少是哪个音调请用傅里叶级数或者变换的方法分析它的谐波分量分别是什么.提示：简单的方法是近似取出一个周期求傅里叶级数但这样明显不准确,由于你应该已经发现基音周期不是整数(这里不允许使用resample函数).复杂些的方法是对整个信号求傅里叶变换(回忆周期性信号 的傅里叶变换),但你可能发现无论你如何提升频域的分辨率,也得不到精确的 包络(应该近似于冲激函数而不是sinc函数),可选的方法是增加时域的数据 量,即再把时域信号重复假设干次,看看这样是否效果好多了请解释之.答：共计221-2=219个周期；基因频率约为8000 / (219 / 9) = 3

17、28.7671Hz. 用FFT寸该信号做DF筏换：10ad( 'Guitar.MAT' );my_wave2proc = repmat(wave2proc, 25, 1);Fs = 8000;Length = 1ength(my_wave2proc);NFFT = 2 " nextpow2(Length);Y = fft(my_wave2proc, NFFT) / Length;amplitude = 2 * abs(Y(1 : NFFT / 2 + 1);frequency = Fs / 2 * linspace(0, 1 ,NFFT / 2 + 1);plot(f

18、requency, amplitude);max_rate_of_grade, max_position = max(amplitude(1 : 100);frequency(max_position)05001000150020002500300035004000算得基频为329.1016Hz,和前面通过周期估算得到的很相近.对照可得,为“调mi.(9)再次载入fmt.wav ,现在要求你写一段程序,自动分析出这段乐曲的音调和节拍！如果你觉得太难就允许手工标定出每个音调的起止时间,再不行你就把每个音调的数据都单独保存成一个文件,然后让MATLAB对这些文件进行批处理.注意：不允许逐一地手工分

19、析音调.编辑音乐文件,推荐使用 CoolEdit编辑软 件.此问参考了学长的版本.看了版本才明白了算法主要分2个步骤：1、分割出单个音符,计算每个音符的时间,2、计算每个音符的频率,把频率转化为音名.分割音符即将一个一个的音符片段从整个曲子中切下.音符的起始都伴随时域上瞬时能量的激增,即时域上波形幅度突然变大.利用这一点特征对音符进行提取.0.01选择小步幅对fmt.wav进行扫描,步幅选为0.01秒最为适宜.具体方法如下：在秒的范围内找极大点,而且是后续能量连续下降的极大点,这个极大点如果满足以下三个条件,那么为有效的音符的起始点：1、这个极大点比前一段的极大点能量高出80%;Pw>1

20、.8*Pwp2、这个极大点比整个乐曲的平均能量高；St>Stp+10003、距前面最近的一个音符的起始点时间差1/8秒以上.Pw>Avg分割的结果：LanHW.S一共分割得30个音,和 CoolEdit的分析根本一致.过度局部Time=diff(StartPoint Len)/8000;EndPoint=StartPoint(2:Num) Len-200;Power=Power/max(Power);STime=Time-200/8000;每个音符的持续时间用于傅立叶变换的区间能量修正(归一)用于傅立叶变换的区间长度计算频率并转化对每个音符内部的一个区间(我选这个音符的起始到结束前

21、的0.025秒,即StartPoint到EndPoint)做快速傅立叶变换,然后在频谱中选出基频.关于自动寻找基频： 一般的,基波不是幅度最大的频率,由于有很多音调的谐波, 甚至非线性谐波的幅度大于基波.一般来说,一次谐波的幅度不太大,不会大于基波幅度的2倍,这样,在有可能出现基波的地方,利用每个点的幅度值除以该点的频率,也就是每个点 对于原点的斜率做为判定基波的标准,经过这样的简化处理, 找到斜率最大的点, 根本上就能够判定基波出现在什么位置.FO=F_fft./omg;%每个点的幅度值除以该点的频率M1 M2=max(FO);%找到斜率最大的点,判定为基波Freq(p)=omg(M2)/2

22、/pi;% 记下频率值Name(p尸Freq2Name(Freq(p);%把频率值转化为音名(关于函数 Freq2Name(Freq)：把频率与十几个标准音的频率作比较,选出最接近的一 个,作为这个频率对应的音名.)display('Score ： ');display(Name);display(Time); %输出结果:#音名时长#音名时长#音名时长#音名时长1-60.16769-40.504417-60.38902510.47182-61.486610-40.74841850.268726-40.53363-70.476411-71.23731940.230527-60.

23、46274-60.45911230.78692030.202128-70.5771520.458413-60.22362120.261029-60.6410630.450014-60.47132210.473130-51.39747-50.466015-60.472823-70.44378-40.44691660.49792420.5484如此,便完成了乐曲的音调和节拍的自动分析.把此结果用函数CreateMusic处理并回放,听到的音乐是很接近原来的fmt.wav的,仅仅有两个左右的音符有偏差.代码如下：fmt=wavread( 'fmt.wav' );%sound(fmt,

24、 8000)Len=length(fmt);Avg=norm(fmt)/sqrt(Len);Num=0;Step=80;Pwp=0;Stp=-1000;St=1;for n=1:Step:Len-StepM=0;St=n;while M=1 & St+Step-1<=LenPw M=max(fmt(St:St+Step-1);St=St+M-1;endif Pw>1.8*Pwp & St>Stp+1000 & Pw>AvgNum=Num+1;StartPoint(Num)=St;Power(Num)=Pw;Stp=St; endPwp=Pw; e

25、ndTime=diff(StartPoint Len)/8000;STime=Time-200/8000;EndPoint=StartPoint(2:Num) Len-200;Power=Power/max(Power);figure;hold on ;plot(fmt);axis(0 length(fmt) 0 0.6);title( ' o ? ? ?);xlabel( 'Time' );ylabel( 'Signal' );for p=1:Numplot(StartPoint(p) StartPoint(p),0 0.6,'r*-'

26、);endFreqLimit=170;for p=1:Numf=fmt(StartPoint(p):EndPoint(p)-1);F_fft,t,omg=FastFourier(f);omg(floor(STime(p)*(4000-FreqLimit):ceil(STime(p)*(4000+FreqLimit)=1e5;FO=F_fft./omg;M1 M2=max(FO);Freq(p)=omg(M2)/2/pi;Name(p)=Freq2Name(Freq(p);enddisplay( 'Score £ o' );display(Name);display(T

27、ime);UnknownMusic=Name;Time;Power;Music=CreateMusic( 'C' ,UnknownMusic, 'Exp' , 'Guitar+' );PlayAndPlot(Music);(10)用(7)计算出来的傅里叶级数再次完成第(4)题,听一听是否像演奏fmt.wav的吉他演奏出来的答：在(7)中计算出来的傅里叶级数中取前四项,后面的谐波并不显著将其略 去.基波和2、3、4次谐波的系数为如下,将其归一后再合成即可.合成的方法 参考了学长的报告.music_score = 5 5 6 2 1 1 6 2;0

28、0 0 0 0 0 -1 0;1 0.5 0.5 2 1 0.5 0.5 2;time = 5;harmonic_coefficient = 0.0538 0.0753 0.0464 0.0482;harmonic_coefficient = harmonic_coefficient /norm(harmonic_coefficient);diationic_scale = 220 * 2 A(3 5 7 8 10 12 14 + 5) / 12);t = 0 : 1 / 8000 : time;music_wave = zeros(1, length(t);time_start = 0;sp

29、eed = 0.5;for n = 1 : length(music_score)frequency = diationic_scale(music_score(1, n) * 2 A music_score(2,n);time_end = time_start + music_score(3, n) * speed ;envelope = GetEnvelope3(time, time_start, time_end);for m = 1 : 4music_wave = music_wave + harmonic_coefficient(m) * sin(2 * pi *m * frequency * t) .* (t >= time_start & t < time_end) .* envelope;endtime_start = time_end;endmusic_wave = music_wave / max(music_wave); sound(music_wave, 8000);(11)也许(9)还不是很像,由于对于一把泛音丰富的吉他而言,不可能每个音调 对应的泛音数量和幅度都相同.但是通过完成第(8)题,你已经提取出fmt.wav中的很多音调,或者说,掌握了每个音调对应的傅里叶级数,大致了解了这把吉他的特征.现在就来演奏一曲?东方红?吧.提