matlab音乐处理合成实验报告

MATLAB高级编程与工程应用

语音合成综合实验

姓名:

班级:

学号:

日期:

1.2.1 简单的合成音乐

(1) 请根据《东方红》片断的简谱和“十二平均律”计算出该片断中各个乐音的频率,在MATLAB 中生成幅度为1 、抽样频率为8kHz 的正弦信号表示这些乐音。请用sound 函数播放每个乐音,听一听音调是否正确。最后用这一系列乐音信号拼出《东方红》片断,注意控制每个乐音持续的时间要符合节拍,用sound 播放你合成的音乐,听起来感觉如何?

由“十二平均律”计算得到各个乐音的频率:

“5”——“C”:523.25Hz

“6”——“D”:587.33Hz

“1”——“F”:349.23Hz

“2”——“G”:392Hz

“6.”频率是“6”的一半:293.66Hz

代码:(project1_1_1.m)

f=8000;

T=1/f;

t8=0:T:1*0.25;

t4=0:T:2*0.25;

t2=0:T:4*0.25;

t1=0:T:8*0.25;

part1=sin(2*pi*523.25*t4);

part2=sin(2*pi*523.25*t8);

part3=sin(2*pi*587.33*t8);

part4=sin(2*pi*392.00*t2);

part5=sin(2*pi*349.23*t4);

part6=sin(2*pi*349.23*t8);

part7=sin(2*pi*293.66*t8);

part8=sin(2*pi*392.00*t2);

total=[part1,part2,part3,part4,part5,part6,part7,part8];

sound(total);

    试听发现,合成后的音乐基本保持了《东方红》的音调,但声音比较沉闷,相邻乐音之间有比较明显的“啪”的杂音。

    (2) 你一定注意到(1) 的乐曲中相邻乐音之间有“啪”的杂声,这是由于相位不连续产生了高频分量。这种噪声严重影响合成音乐的质量,丧失真实感。为了消除它,我们可以用图1.5 所示包络修正每个乐音,以保证在乐音的邻接处信号幅度为零。此外建议用指数衰减的包络来表示。

    首先尝试用折线包络,编写函数envelope_line.m生成所需折线:

function envelope = envelope_line(t)

envelope(1:floor(t/8)) = linspace(0,1,floor(t/8));

envelope(floor(t/8)+1:floor(t/4)) = linspace(1,0.5,floor(t/4)-floor(t/8)));

envelope(floor(t/4)+1:floor(3*t/4)) = linspace(0.5,0.5,floor(t*3/4)-floor(t/4));

envelope(floor(3*t/4)+1:t) = linspace(0.5,0,floor(t)-floor(3*t/4));

对project1_1_1.m中的部分代码进行修改,调用envelope_line实现折线包络:(project1_1_2.m)

part1=sin(2*pi*523.25*t4).*envelope_line(t4);

part2=sin(2*pi*523.25*t8).*envelope_line(t8);

part3=sin(2*pi*587.33*t8).*envelope_line(t8);

part4=sin(2*pi*392.00*t2).*envelope_line(t2);

part5=sin(2*pi*349.23*t4).*envelope_line(t4);

part6=sin(2*pi*349.23*t8).*envelope_line(t8);

part7=sin(2*pi*293.66*t8).*envelope_line(t8);

part8=sin(2*pi*392.00*t2).*envelope_line(t2);

    试听结果杂音明显减少,但声音还比较生硬,尝试用指数衰减的包络来表示。也是对project1_1_1.m中的部分代码进行修改(project1_1_2.m)

part1=sin(2*pi*523.25*t4).*exp(-2*t4);

part2=sin(2*pi*523.25*t8).*exp(4*t8);

part3=sin(2*pi*587.33*t8).*exp(4*t8);

part4=sin(2*pi*392.00*t2).*exp(1*t2);

part5=sin(2*pi*349.23*t4).*exp(2*t4);

part6=sin(2*pi*349.23*t8).*exp(4*t8);

part7=sin(2*pi*293.66*t8).*exp(4*t8);

part8=sin(2*pi*392.00*t2).*exp(1*t2);

    试听结果比较好,声音很圆润,但仔细听的话也能发现杂音。

    (3) 请用最简单的方法将(2) 中的音乐分别升高和降低一个八度。(提示:音乐播放的时间可以变化)再难一些,请用resample 函数(也可以用interp 和decimate 函数)将上述音乐升高半个音阶。(提示:视计算复杂度,不必特别精确)

方法一:

在正弦信号内添加系数改变其频率。

升高八度:

part1=sin(4*pi*523.25*t4).*exp(-4*t4);

part2=sin(4*pi*523.25*t8).*exp(8*t8);

part3=sin(4*pi*587.33*t8).*exp(8*t8);

part4=sin(4*pi*392.00*t2).*exp(2*t2);

part5=sin(4*pi*349.23*t4).*exp(4*t4);

part6=sin(4*pi*349.23*t8).*exp(8*t8);

part7=sin(4*pi*293.66*t8).*exp(8*t8);

part8=sin(4*pi*392.00*t2).*exp(2*t2);

high8_total=[part1,part2,part3,part4,part5,part6,part7,part8];

sound(high8_total);

降低八度:

part1=sin(1*pi*523.25*t4).*exp(-1*t4);

part2=sin(1*pi*523.25*t8).*exp(2*t8);

part3=sin(1*pi*587.33*t8).*exp(2*t8);

part4=sin(1*pi*392.00*t2).*exp(0.5*t2);

part5=sin(1*pi*349.23*t4).*exp(1*t4);

part6=sin(1*pi*349.23*t8).*exp(2*t8);

part7=sin(1*pi*293.66*t8).*exp(2*t8);

part8=sin(1*pi*392.00*t2).*exp(0.5*t2);

low8_total=[part1,part2,part3,part4,part5,part6,part7,part8];

sound(low8_total);

方法二:

       直接调用resample函数。

升高八度:

high8_total=resample(total,1,2);

降低八度:

low8_total=resample(total,2,1);

升高半个音阶:

根据“十二平均律”中2^(1/12)=1.06

highhalf_total=resample(total,1,2);

    (4) 试着在(2) 的音乐中增加一些谐波分量,听一听音乐是否更有“厚度”了?注意谐波分量的能量要小,否则掩盖住基音反而听不清音调了。(如果选择基波幅度为1 ,二次谐波幅度0:2 ,三次谐波幅度0:3 ,听起来像不像象风琴?)

(project1_1_4.m)

part1=sin(2*pi*523.25*t4).*exp(-2*t4)+0.2*sin(4*pi*523.25*t4).*exp(-2*t4)+0.3*sin(6*pi*523.25*t4).*exp(-2*t4);

part2=sin(2*pi*523.25*t8).*exp(4*t8)+0.2*sin(4*pi*523.25*t8).*exp(4*t8)+0.3*sin(6*pi*523.25*t8).*exp(4*t8);

part3=sin(2*pi*587.33*t8).*exp(4*t8)+0.2*sin(4*pi*587.33*t8).*exp(4*t8)+0.3*sin(6*pi*587.33*t8).*exp(4*t8);

part4=sin(2*pi*392.00*t2).*exp(1*t2)+0.2*sin(4*pi*392.00*t2).*exp(1*t2)+0.3*sin(6*pi*392.00*t2).*exp(1*t2);

part5=sin(2*pi*349.23*t4).*exp(2*t4)+0.2*sin(4*pi*349.23*t4).*exp(2*t4)+0.3*sin(6*pi*349.23*t4).*exp(2*t4);

part6=sin(2*pi*349.23*t8).*exp(4*t8)+0.2*sin(4*pi*349.23*t8).*exp(4*t8)+0.3*sin(6*pi*349.23*t8).*exp(4*t8);

part7=sin(2*pi*293.66*t8).*exp(4*t8)+0.2*sin(4*pi*293.66*t8).*exp(4*t8)+0.3*sin(6*pi*293.66*t8).*exp(4*t8);

part8=sin(2*pi*392.00*t2).*exp(1*t2)+0.2*sin(4*pi*392.00*t2).*exp(1*t2)+0.3*sin(6*pi*392.00*t2).*exp(1*t2);

试听结果确实变得醇厚了,有风琴的感觉。

    (5) 自选其它音乐合成,例如贝多芬第五交响乐的开头两小节。

    我选取了一首儿歌《粉刷匠》(project1_1_5.m)

粉刷匠(波兰儿歌)_简谱_搜谱网.jpg

f=8000;

T=1/f;

t8=0:T:1*0.25;

t4=0:T:2*0.25;

t2=0:T:4*0.25;

t1=0:T:8*0.25;

part1=sin(2*pi*523.25*t8).*exp(4*t8);

part2=sin(2*pi*440.00*t8).*exp(4*t8);

part3=sin(2*pi*523.25*t8).*exp(4*t8);

part4=sin(2*pi*440.00*t8).*exp(4*t8);

part5=sin(2*pi*523.25*t8).*exp(4*t8);

part6=sin(2*pi*440.00*t8).*exp(4*t8);

part7=sin(2*pi*349.23*t4).*exp(2*t4);

part8=sin(2*pi*392.00*t8).*exp(4*t8);

part9=sin(2*pi*493.88*t8).*exp(4*t8);

part10=sin(2*pi*440.00*t8).*exp(4*t8);

part11=sin(2*pi*392.00*t8).*exp(4*t8);

part12=sin(2*pi*523.25*t2).*exp(1*t2);

total=[part1,part2,part3,part4,part5,part6,part7,part8,part9,part10,part11,part12];

plot(total);

sound(total);

1.2.2 用傅里叶级数分析音乐

(6) 先用wavread 函数载入光盘中的fmt.wav 文件,播放出来听听效果如何?是否比刚才的合成音乐真实多了?

(project1_1_6.m)

x=wavread('fmt.wav');

plot(x);

sound(x);

试听结果的确真实许多。

(7) 你知道待处理的wave2proc 是如何从真实值realwave 中得到的么?这个预处理过程可以去除真实乐曲中的非线性谐波和噪声,对于正确分析音调是非常重要的。提示:从时域做,可以继续使用resample 函数。

    这题开始时没有理解题目要求,不知道如何入手,于是在CSDN上查找帮助。要去除非线性杂波和噪声,首先要找出它们与真实乐曲的区别。非线性杂波和噪声都是随机产生的,不具有周期性,因此要在周期性的乐音中将其去除,可以考虑将真实乐曲多次叠加再取平均值。乐音可以想象,叠加次数越多,最后得到的平均后的乐音越具有周期性。

    读出realwave波形,发现采样点为243,重复10次,可得周期为24.3,所以可以将其延长至10倍,这样,周期就是整数。

(project1_1_7.m)

load('guitar.mat');

l=length(realwave);

add=resample(realwave,10,1);

average=zeros(1,l);

for m=1:10

average=average+(add((m-1)*l+1:m*l))';

end

add=[average,average,average,average,average,average,average,average,average,average];

average=resample(add/10,1,10);

figure;

subplot(3,1,1);

plot(wave2proc);

title('wave2proc' );

subplot(3,1,2);

plot(average);

title('average');

subplot(3,1,3);

plot(average-wave2proc');

title('average-wave2proc');

由图中前两个波形可以看出,所得average 波形与wave2proc波形几乎重合。但由第三个波形看出两者还是有细微差别的。

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

开始我的代码是这样的:

F=8000;

l1=length(wave2proc);

l=10*l1;

[t omg FT IFT]=prefourier([0,(l-1)/F],l,[0,20000],10000);

r=[wave2proc;wave2proc;wave2proc;wave2proc;wave2proc;wave2proc;wave2proc;wave2proc;wave2proc;wave2proc];

R=FT*r;

R1=abs(R);

plot(omg,R1);

set(gca,'XGrid','on','YGrid','on');

但是对于prefourier函数的调用总是出错,(l-1)/F这一项不能为double。估计是MATLAB版本问题,于是换了一种fourier变换方法。

先将信号重复100遍,然后对整个信号求傅里叶变换,然后作图即可求得x的基频。

(project1_1_8.m)

load('guitar.mat');

Times=100;

s=repmat(wave2proc,Times,1);

L=243*Times;

N=L;

T=L/8000;

OMG=N/T*2*pi;

t=linspace(-T/2,T/2-T/L,L)';

omg=linspace(-OMG/2,OMG/2-OMG/L,L)';

f1=s.*exp(-j*(-OMG/2)*t);

F1=T*exp(j*(-OMG/2)*(-T/2))/N*fft(f1);

plot(omg,abs(F1),'b-');

如图所示: 基频频率为:329.2Hz,幅值为0.08264;

二次谐波:幅值为0.1204;

三次谐波:幅值为0.007923; 四次谐波:幅值为0.0909;

我认为重复的好处是第一采样的点数多了,这样得到的信息量会增多。另外重复周期多了,就越接近周期信号,这样的话,傅里叶变换得到的图会比较清晰,像冲击函数。而如果只是对wave2proc进行傅里叶变换的话就会看到每一个谐波分量附近都有很明显是sa或者升余弦这样的形状。

(9) 再次载入fmt.wav ,现在要求你写一段程序,自动分析出这段乐曲的音调和节拍!如果你觉得太难就允许手工标定出每个音调的起止时间,再不行你就把每个音调的数据都单独保存成一个文件,然后让MATLAB 对这些文件进行批理。注意:不允许逐一地手工分析音调。编辑音乐文件,推荐使用\CoolEdit" 编辑软件。

总体思路:先将每个音分开来,然后对同音调的这段音乐作傅里叶变换,然后找到基频。(同时为了后面方便,找到谐波的系数。)

    在对音乐分段的时候,我利用的是音乐包络的突变性。通过max函数对波形取包络,这样可以比较明显的看出那边的变化率比较大。对变化率设定一个范围(在这里引入了不少误差)借此给音频信号分段。分出来的段落和人耳辨别得到的基本一致,但是有若干吉他的颤音被认为是弹奏了。这样的误差如果对音调做一点处理应当是能够消除的。(向同学请教了分段的根据)

在找基频的时候,采用傅里叶变换和上题类似。一开始比较简单的认为最大值即是基频,因为在6.2.1-4的题目中看到谐波分量比较小这一说法。后来参考了前面的题的傅里叶变换,发现谐波分量可以大于基频分量。因此程序做了修改,因为时间紧张,所以稍显仓促。大致的想法还是利用了max函数。因为最大值一定是谐波分量或者基频分量。考虑到这么大分量的谐波分量不会很大,一般都是2-3次谐波,最多是4次谐波,所以采用试探的方法。基频分量不可能小到无法和谐波分量相比,所以对基频分量和最大分量的比例设下限来寻找基频。当找到一个大于下限比的值时,如果它和最大分量的比值几乎为整数,则可认为是基频。

(project1_1_9.m)

y=wavread('fmt.wav');

n=length(y);

k=100;

a1=floor(n/k);

x1=zeros(a1+1,3);

x2=zeros(a1+1,8);

m=1;

x2(1,1)=0;

for a=0:a1

    if(a==a1)

        [x1(a+1,1),x1(a+1,2)]=max(y(k*a+1:n));

    else

        [x1(a+1,1),x1(a+1,2)]=max(y(k*a+1:k*(a+1)));

    end

end

for a=1:a1

    if(x1(a+1,1)-x1(a,1)>0.06)

        x1(a+1,3)=1;

        x2(m,2)=((a+1)*k+x1(a+1,2))/8000;

        m=m+1;

        x2(m,1)=x2(m-1,2);

    end

end

for a=1:a1

    if(x2(a,2)~=0)

        k1=floor(x2(a,1)*8000)+1;

        k2=floor(x2(a,2)*8000);

        s=repmat(y(k1:k2),100,1);

        L=length(s);

        N=L;

        T=L/8000;

        OMG=N/T*2*pi;

        t=linspace(-T/2,T/2-T/L,L)';

        omg=linspace(-OMG/2,OMG/2-OMG/L,L)';

        f1=s.*exp(-j*(-OMG/2)*t);

        F=T*exp(j*(-OMG/2)*(-T/2))/N*fft(f1);

        F1=abs(F);

        [C,I]=max(F1(N/2:N));

        for A=1:I

            if(F1(A+N/2)>C/3)

                if(abs(I/A-5)<0.01)

                    x2(a,3)=A*8000/L;

                    x2(a,4)=F1(A+N/2);

                    x2(a,5)=F1(2*A+N/2);

                    x2(a,6)=F1(3*A+N/2);

                    x2(a,7)=F1(4*A+N/2);

                    x2(a,8)=C;

                    break;

                else if(abs(I/A-4)<0.01)

                        x2(a,3)=A*8000/L;

                        x2(a,4)=F1(A+N/2);

                        x2(a,5)=F1(2*A+N/2);

                        x2(a,6)=F1(3*A+N/2);

                        x2(a,7)=C;

                        x2(a,8)=F1(5*A+N/2);

                        break;

                    else if(abs(I/A-3)<0.01)

                            x2(a,3)=A*8000/L;

                            x2(a,4)=F1(A+N/2);

                            x2(a,5)=F1(2*A+N/2);

                            x2(a,6)=C;

                            x2(a,7)=F1(4*A+N/2);

                            x2(a,8)=F1(5*A+N/2);

                            break;

                        else if(abs(I/A-2)<0.01)

                                x2(a,3)=A*8000/L;

                                x2(a,4)=F1(A+N/2);

                                x2(a,5)=C;

                                x2(a,6)=F1(3*A+N/2);

                                x2(a,7)=F1(4*A+N/2);

                                x2(a,8)=F1(5*A+N/2);

                                break;

                            end

                        end

                    end

                end

            end

        end

        if(A==I)

            x2(a,3)=I*8000/L;

            x2(a,4)=C;

            x2(a,5)=F1(2*I+N/2);

            x2(a,6)=F1(3*I+N/2);

            x2(a,7)=F1(4*I+N/2);

            x2(a,8)=F1(5*I+N/2);

        end

    end

end

1.2.3 基于傅里叶级数的合成音乐

现在进入了合成音乐的高级境界,我们要用演奏fmt.wav 的吉他合成出一段《东方红》。

(10) 用(7) 计算出来的傅里叶级数再次完成第(4) 题,听一听是否像演奏fmt.wav 的吉他演奏出来的?

(project1_1_10.m)

fz=8000;

t=[0:1/fz:4-1/fz];

t1=t(1:0.5*fz);

t2=t(0.5*fz+1:0.75*fz);

t3=t(0.75*fz+1:1*fz);

t4=t(1*fz+1:2*fz);

t5=t(2*fz+1:2.5*fz);

t6=t(2.5*fz+1:2.75*fz);

t7=t(2.75*fz+1:3*fz);

t8=t(3*fz+1:4*fz);

p1=523.25;

p2=523.25;

p3=587.33;

p4=392;

p5=349.23;

p6=349.23;

p7=293.66;

p8=392;

part1=exp(-2*t1).*(0.413*sin(2*pi*p1*t1)+0.6016*sin(4*pi*p1*t1)+0.3957*sin(6*pi*p1*t1)+0.4539*sin(8*pi*p1*t1));

part2=exp(4*(0.5-t2)).*(0.413*sin(2*pi*p2*t2)+0.6016*sin(4*pi*p2*t2)+0.3957*sin(6*pi*p2*t2)+0.4539*sin(8*pi*p2*t2));

part3=exp(4*(0.75-t3)).*(0.413*sin(2*pi*p3*t3)+0.6016*sin(4*pi*p3*t3)+0.3957*sin(6*pi*p3*t3)+0.4539*sin(8*pi*p3*t3));

part4=exp(1*(1-t4)).*(0.413*sin(2*pi*p4*t4)+0.6016*sin(4*pi*p4*t4)+0.3957*sin(6*pi*p4*t4)+0.4539*sin(8*pi*p4*t4));

part5=exp(2*(2-t5)).*(0.413*sin(2*pi*p5*t5)+0.6016*sin(4*pi*p5*t5)+0.3957*sin(6*pi*p5*t5)+0.4539*sin(8*pi*p5*t5));

part6=exp(4*(2.5-t6)).*(0.413*sin(2*pi*p6*t6)+0.6016*sin(4*pi*p6*t6)+0.3957*sin(6*pi*p6*t6)+0.4539*sin(8*pi*p6*t6));

part7=exp(4*(2.75-t7)).*(0.413*sin(2*pi*p7*t7)+0.6016*sin(4*pi*p7*t7)+0.3957*sin(6*pi*p7*t7)+0.4539*sin(8*pi*p7*t7));

part8=exp(1*(3-t8)).*(0.413*sin(2*pi*p8*t8)+0.6016*sin(4*pi*p8*t8)+0.3957*sin(6*pi*p8*t8)+0.4539*sin(8*pi*p8*t8));

total=[part1,part2,part3,part4,part5,part6,part7,part8];

sound(total);

    我发誓我是想弹一段吉他出来,但试听的结果像是钢琴,不过还是挺好听的,至少是个乐器。

(11) 也许(9) 还不是很像,因为对于一把泛音丰富的吉他而言,不可能每个音调对应的泛音数量和幅度都相同。但是通过完成第(8) 题,你已经提取出fmt.wav 中的很多音调,或者说,掌握了每个音调对应的傅里叶级数,大致了解了这把吉他的特征。现在就来演奏一曲《东方红》吧。提示:如果还是音调信息不够,那就利用相邻音调的信息近似好了,毕竟可以假设吉他的频响是连续变化的。

(project1_1_11.m)

fz=8000;

t=[0:1/fz:4-1/fz];

t1=t(1:0.5*fz);

t2=t(0.5*fz+1:0.75*fz);

t3=t(0.75*fz+1:1*fz);

t4=t(1*fz+1:2*fz);

t5=t(2*fz+1:2.5*fz);

t6=t(2.5*fz+1:2.75*fz);

t7=t(2.75*fz+1:3*fz);

t8=t(3*fz+1:4*fz);

p1=523.25;

p2=523.25;

p3=587.33;

p4=392;

p5=349.23;

p6=349.23;

p7=293.66;

p8=392;

part1=exp(-2*t1).*(sin(2*pi*p1*t1)+0.73339*sin(4*pi*p1*t1));

part2=exp(4*(0.5-t2)).*(sin(2*pi*p2*t2)+0.73339*sin(4*pi*p2*t2));

part3=exp(4*(0.75-t3)).*(sin(2*pi*p3*t3)+0.33841*sin(6*pi*p3*t3));

part4=exp(1*(1-t4)).*(sin(2*pi*p4*t4)+0.6842*sin(4*pi*p4*t4));

part5=exp(2*(2-t5)).*(0.5*sin(2*pi*p5*t5)+0.66482*sin(4*pi*p5*t5)+0.4840*sin(6*pi*p5*t5)+0.2654*sin(8*pi*p5*t5));

part6=exp(4*(2.5-t6)).*(0.5*sin(2*pi*p6*t6)+0.66482*sin(4*pi*p6*t6)+0.4840*sin(6*pi*p6*t6)+0.2654*sin(8*pi*p6*t6));

part7=exp(4*(2.75-t7)).*(sin(2*pi*p7*t7)+0.96774*sin(4*pi*p7*t7));

part8=exp(1*(3-t8)).*(sin(2*pi*p8*t8)+0.6842*sin(4*pi*p8*t8));

total=[part1,part2,part3,part4,part5,part6,part7,part8];

sound(total);

    这次是只有一根弦的吉他?有吉他的感觉,但声音不够丰满。我想这是因为泛音需要很多频率的叠加,而不仅仅是几个谐波加一个基波。

实验收获及感想:

         这次大作业总体上比语音合成要难,花费的时间也多。其中很大一笔时间用在纠错上,最常见的错误是函数的参数不合要求。有时明明按例题给参数,也会报错,仔细阅读函数说明会发现不同版本MATLAB对同一函数的描述也是不一样的,甚至对参数的要求也不同。

         对于一段比较长的代码要分开编写,及时检测,勤代码中插入plot语句,及时发现错误并纠正。对于不同方法编辑出的同一段乐音,要多对比试听,发现不同因素对乐音的影响。

         这次大作业花费的时间很长,但使我真正掌握了傅里叶变换的应用,还是很有收获的。

相关推荐