信号与系统
- 综合实验之音乐合成
(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*omg5*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);
听的时候发现在相邻乐音之间有杂音,这是由于相位不连续造成的。
(2) 你一定注意到(1) 的乐曲中相邻乐音之间有“啪”的杂声,这是由于相位不连续产生了高频分量。这种噪声严重影响合成音乐的质量,丧失真实感。为了消除它,我们可以用图1.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*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];
sound(m);
第一次我采用的指数衰减没有时间前面的系数,即每个都只乘exp(t),没有系数;后来根据不同节拍,更改了不同的衰减系数,这样声音听起来感觉更加圆润。
(3) 请用最简单的
将(2) 中的音乐分别升高和降低一个八度。(提示:音乐播放的时间可以变化)再难一些,请用resample 函数(也可以用interp 和decimate 函数)将上述音乐升高半个音阶。(提示:视计算复杂度,不必特别精确)
答:最简单的方法是直接更改抽样频率f。将f从8K改为4K,则升高一个八度,并且播放速度增快了一倍;将f从8k改为16k,则降低一个八度,速度也变慢了一倍。
升高半个音阶,只须在(2)代码最后加一句resample(m,1000,1059)即可。
f =8000; %改为4000或者16000
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*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,1000,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*omg5*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*t2).*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*omg6l*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;
t2=[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*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(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 中得到的么?这个预处理过程可以去除真实乐曲中的非线性谐波和噪声,对于正确分析音调是非常重要的。提示:从时域做,可以继续使用resample 函数。
realwave中的波形有十个周期,要除去其中的噪声可以采用时域求均值的方法。步骤如下:首先用resample函数将其采样率增大为十倍,再等分十份取平均,然后重复这个平均后的波形十次,还原其长度,最后在用resample函数还原采样率到1/10。
代码如下:
load('Guitar.MAT');
wavetemp = zeros(length(realwave), 1);
waveresampled = resample(realwave,10,1);
for n = 1 : 10
wavetemp = wavetemp + waveresampled((n - 1) * length(realwave) + 1 : n * length(realwave)) / 10;
end
mywave2proc = 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);
(8) 这段音乐的基频是多少?是哪个音调?请用傅里叶级数或者变换的方法分析它的谐波分量分别是什么。提示:简单的方法是近似取出一个周期求傅里叶级数但这样明显不准确,因为你应该已经发现基音周期不是整数(这里不允许使用resample 函数)。复杂些的方法是对整个信号求傅里叶变换(回忆周期性信号的傅里叶变换),但你可能发现无论你如何提高频域的分辨率,也得不到精确的包络(应该近似于冲激函数而不是sinc 函数),可选的方法是增加时域的数据量,即再把时域信号重复若干次,看看这样是否效果好多了?请解释之。
答:共计221-2=219个周期;基因频率约为8000 / (219 / 9) = 328.7671Hz.
用FFT对该信号做DFT变换:
load('Guitar.MAT');
my_wave2proc = repmat(wave2proc, 25, 1);
Fs = 8000;
Length = length(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(frequency, amplitude);
[max_rate_of_grade, max_position] = max(amplitude(1 : 100));
frequency(max_position)
算得基频为329.1016Hz,和前面通过周期估算得到的很相近。
对照可得,为C大调mi。
(9) 再次载入fmt.wav ,现在要求你写一段程序,自动分析出这段乐曲的音调和节拍!如果你觉得太难就允许手工标定出每个音调的起止时间,再不行你就把每个音调的数据都单独保存成一个文件,然后让MATLAB 对这些文件进行批处理。注意:不允许逐一地手工分析音调。编辑音乐文件,推荐使用\CoolEdit" 编辑软件。
此问参考了学长的版本。看了版本才明白了算法。
主要分2个步骤:
1、分割出单个音符,计算每个音符的时间,2、计算每个音符的频率,把频率转化为音名。
分割音符
即将一个一个的音符片段从整个曲子中切下。音符的起始都伴随时域上瞬时能量的激增,即时域上波形幅度突然变大。利用这一点特征对音符进行提取。
选择小步幅对fmt.wav进行扫描,步幅选为0.01秒最为适宜。具体方法如下:在0.01秒的范围内找极大点,而且是后续能量连续下降的极大点,这个极大点如果满足以下三个条件,则为有效的音符的起始点:
1、 这个极大点比前一段的极大点能量高出80%;
Pw>1.8*Pwp
2、 这个极大点比整个乐曲的平均能量高;
St>Stp+1000
3、 距前面最近的一个音符的起始点时间差1/8秒以上。
Pw>Avg
分割的结果:
一共分割得30个音,和CoolEdit的分析基本一致。
过度部分
每个音符的持续时间
Time=diff([StartPoint Len])/8000;
用于傅立叶变换的区间
EndPoint=[StartPoint(2:Num) Len]-200;
能量修正(归一)
Power=Power/max(Power);
用于傅立叶变换的区间长度
STime=Time-200/8000;
计算频率并转化
对每个音符内部的一个区间(我选这个音符的起始到结束前的0.025秒,即StartPoint到EndPoint)做快速傅立叶变换,然后在频谱中选出基频。
关于自动寻找基频:一般的,基波不是幅度最大的频率,因为有很多音调的谐波,甚至非线性谐波的幅度大于基波。一般来说,一次谐波的幅度不太大,不会大于基波幅度的2倍,这样,在有可能出现基波的地方,利用每个点的幅度值除以该点的频率,也就是每个点对于原点的斜率做为判定基波的
,经过这样的简化处理,找到斜率最大的点,基本上就能够判定基波出现在什么位置。
FO=F_fft./omg;
%每个点的幅度值除以该点的频率
[M1 M2]=max(FO);
%找到斜率最大的点,判定为基波
Freq(p)=omg(M2)/2/pi;
%记下频率值
Name(p)=Freq2Name(Freq(p));
%把频率值转化为音名
(关于函数Freq2Name(Freq):把频率与十几个标准音的频率作比较,选出最接近的一个,作为这个频率对应的音名。)
display('Score:');display(Name);display(Time);
%输出结果:
#
音名
时长
#
音名
时长
#
音名
时长
#
音名
时长
1
-6
0.1676
9
-4
0.5044
17
-6
0.3890
25
1
0.4718
2
-6
1.4866
10
-4
0.7484
18
5
0.2687
26
-4
0.5336
3
-7
0.4764
11
-7
1.2373
19
4
0.2305
27
-6
0.4627
4
-6
0.4591
12
3
0.7869
20
3
0.2010
28
-7
0.5771
5
2
0.4584
13
-6
0.2236
21
2
0.2610
29
-6
0.6410
6
3
0.4500
14
-6
0.4713
22
1
0.4731
30
-5
1.3974
7
-5
0.4660
15
-6
0.4728
23
-7
0.4437
8
-4
0.4469
16
6
0.4979
24
2
0.5484
如此,便完成了乐曲的音调和节拍的自动分析。把此结果用函数CreateMusic处理并回放,听到的音乐是很接近原来的fmt.wav的,仅仅有两个左右的音符有偏差。
代码如下:
fmt=wavread('fmt.wav');
%sound(fmt, 8000)
Len=length(fmt);
Avg=norm(fmt)/sqrt(Len);
Num=0;Step=80;Pwp=0;Stp=-1000;St=1;
for n=1:Step:Len-Step
M=0;St=n;
while M~=1 & St+Step-1<=Len
[Pw M]=max(fmt(St:St+Step-1));
St=St+M-1;
end
if Pw>1.8*Pwp & St>Stp+1000 & Pw>Avg
Num=Num+1;
StartPoint(Num)=St;
Power(Num)=Pw;
Stp=St;
end
Pwp=Pw;
end
Time=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('Òô·û·Ö¸î');
xlabel('Time');ylabel('Signal');
for p=1:Num
plot([StartPoint(p) StartPoint(p)],[0 0.6],'r*-') ;
end
FreqLimit=170;
for p=1:Num
f=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));
end
display('Score£º');display(Name);display(Time);
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 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.^(([3 5 7 8 10 12 14] + 5) / 12)];
t = [0 : 1 / 8000 : time];
music_wave = zeros(1, length(t));
time_start = 0;
speed = 0.5;
for n = 1 : length(music_score)
frequency = diationic_scale(music_score(1, n)) * 2 ^ music_score(2, n);
time_end = time_start + music_score(3, n) * speed ;
envelope = GetEnvelope3(time, time_start, time_end);
for m = 1 : 4
music_wave = music_wave + harmonic_coefficient(m) * sin(2 * pi * m * frequency * t) .* (t >= time_start & t < time_end) .* envelope;
end
time_start = time_end;
end
music_wave = music_wave / max(music_wave);
sound(music_wave, 8000);
(11) 也许(9) 还不是很像,因为对于一把泛音丰富的吉他而言,不可能每个音调对应的泛音数量和幅度都相同。但是通过完成第(8) 题,你已经提取出fmt.wav 中的很多音调,或者说,掌握了每个音调对应的傅里叶级数,大致了解了这把吉他的特征。现在就来演奏一曲《东方红》吧。提示:如果还是音调信息不够,那就利用相邻音调的信息近似好了,毕竟可以假设吉他的频响是连续变化的。
答:和上题类似。
music_score = [5 5 6 2 1 1 6 2;
0 0 0 0 0 0 -1 0;
1 0.5 0.5 2 1 0.5 0.5 2];
time = 5;
load('HarmonicCoefficients.mat');
diationic_scale = [220 * 2.^(([3 5 7 8 10 12 14] + 5) / 12)];
t = [0 : 1 / 8000 : time];
music_wave = zeros(1, length(t));
time_start = 0;
speed = 0.5;
for n = 1 : length(music_score)
frequency = diationic_scale(music_score(1, n)) * 2 ^ music_score(2, n);
time_end = time_start + music_score(3, n) * speed ;
envelope = GetEnvelope3(time, time_start, time_end);
musical_name_index = round(12 * log(frequency / 220) / log(2)) + 13;
for m = 1 : 7
music_wave = music_wave + harmonic_coefficient(musical_name_index, m) * sin(2 * pi * m * frequency * t) .* (t >= time_start & t < time_end) .* envelope;
end
time_start = time_end;
end
music_wave = music_wave / max(music_wave);plot(music_wave);
sound(music_wave, 8000);
两种吉他版本听起来都更加清脆悦耳一些。
实验
:
这次试验总体感觉要比第一次难写。一开始我不太会对于音乐信号进行傅里叶变化并且进行后续分析。后来才考了学长的报告,以及和同学讨论后才懂得对于这些函数的运用方法。
最后由于时间以及能力所限,并没有做图形界面。