基于MATLAB的音乐合成和处理

年夜学结业课程设计

工具/原料

MATLAB
MATLAB2016A

方式/步骤

1
基于MATLAB的音乐旌旗灯号合当作与处置
摘   要

本设计共有三部门：1.简单的音乐合当作；2.用傅里叶变换阐发音乐；3.基于傅里叶级数的音乐合当作。本设计采用MATLAB软件仿真来实现。首先,经由过程对音乐旌旗灯号的采样、抽取、调制、解调等多种处置过程的理论阐发和MATLAB求得这段音乐的基频、谐波分量、等数据;然后,经由过程对乐理的研究,按照阐发中求得的数据编写法式,进行基于傅里叶阐发的音乐合当作设计，并设计了图形用户界面；最后加强软件编程实现能力息争决现实问题的能力。
1简单的合当作音乐
1.1 乐理常识介绍
　　
    乐音的根基特征可以用基波频率、谐波频谱和包络波形3个方面来描述。
基波频率：每个指定调子的唱名都对应固定的基波旌旗灯号频率。所谓唱名是指常日读曲谱唱出的1(do)、2(re)、3(mi)… … ,每个唱名并未固定基波频率。当指定乐曲的调子时才知道此时唱名对应的频率值。如C调“ 1”的基波频率为261.63HZ,F调“1”的基波频率为349.23HZ,F调“ 5”的基波频率为523.25HZ。
谐波频谱:在这七个音符中有一个纪律，就是3（mi）到4（fa），7（si）到高音1（do）是半音。在吉他上是相邻的两个品为半音，好比一弦1品是3（mi），那么一弦2品就是4（fa）；在吉他上隔一品是全音，好比一弦1品是1（do），那么一弦3品就是2（re），中心隔了1品。包络波形:分歧类型的乐器,包络外形也不不异。在音乐合当作尝试中,为简化编程描述,凡是把复杂的包络函数用少量直线近似。于是,乐音波形的包络呈拆线。有时为了包管在乐音的邻接处旌旗灯号幅度为零,也可以用指数衰减的包络来暗示,这也是最简单的法子。
1.2 操纵MATLAB实现音乐合当作
本设计采用扬基杜德尔小曲作
按照《扬基杜德尔》第一末节的简谱和十二平均律计较出该末节每个乐音的频率，在MATLAB中生当作幅度为1，抽样频率为8kHz的正弦旌旗灯号暗示这些乐音，用sound播放合当作的音乐。而在MATLAB中暗示乐音所用的抽样频率为fs=8000Hz，也就是所1s钟内有8000个点，抽样点数的几多就可暗示出每个乐音的持续时候的长短。用一个行标的目的量来存储这段音乐对应的抽样点，在用sound函数播放即可。
下为在MATLAB中编写法式
clc
clear
fs=44100;
t=0:1/fs:0.5;
c3_2=key(48, 2, fs);
d3_2=key(50, 2, fs);
e3_2=key(52, 2, fs);
f3_2=key(53, 2, fs);
g3_2=key(55, 2, fs);
a3_2=key(57, 2, fs);
b3_2=key(59, 2, fs);
c3_4=key(48, 4, fs);
d3_4=key(50,4, fs);
e3_4=key(52, 4, fs);
f3_4=key(53, 4, fs);
g3_4=key(55, 4, fs);
a3_4=key(57, 4, fs);
b3_4=key(59, 4, fs);
c3_8=key(48, 8, fs);
d3_8=key(50,8, fs);
e3_8=key(52, 8, fs);
f3_8=key(53, 8, fs);
g3_8=key(55, 8, fs);
a3_8=key(57, 8, fs);
b3_8=key(59, 8, fs);
c3_16=key(48, 16, fs);
d3_16=key(50,16, fs);
e3_16=key(52, 16, fs);
f3_16=key(53, 16, fs);
g3_16=key(55,16, fs);
a3_16=key(57, 16, fs);
b3_16=key(59, 16, fs);
c4_2=key(60, 2, fs);
d4_2=key(62, 2, fs);
e4_2=key(64, 2, fs);
f4_2=key(65, 2, fs);
g4_2=key(67, 2, fs);
a4_2=key(69, 2, fs);
b4_2=key(71, 2, fs);
c4_4=key(60, 4, fs);
d4_4=key(62, 4, fs);
e4_4=key(64, 4, fs);
f4_4=key(65, 4, fs);
g4_4=key(67, 4, fs);
a4_4=key(69, 4, fs);
b4_4=key(71, 4, fs);
c4_8=key(60, 8, fs);
d4_8=key(62, 8, fs);
e4_8=key(64, 8, fs);
f4_8=key(65, 8, fs);
g4_8=key(67, 8, fs);
a4_8=key(69, 8, fs);
b4_8=key(71, 8, fs);
c4_16=key(60, 16, fs);
d4_16=key(62,16, fs);
e4_16=key(64, 16, fs);
f4_16=key(65, 16, fs);
g4_16=key(67, 16, fs);
a4_16=key(69, 16, fs);
b4_16=key(71, 16, fs);
part1=[g3_4 c4_8 c4_8 d4_8 e4_8 c4_8 e4_8 d4_8 b3_8];
part2=[c4_8 c4_8 d4_8 e4_8 c4_4 b3_8 g3_8 ];
part3=[c4_8 c4_8 d4_8 e4_8 f4_8 e4_8 d4_8 c4_8 b3_8 g3_8 a3_8 b3_8 c4_4 c4_4];
part4=[a3_8 b3_16 a3_8 g3_8 a3_8 b3_8 c4_4];
part5=[g3_8 a3_16 g3_8 f3_8 e3_8 f3_8 g3_4];
part6=[a3_8 b3_16 a3_8 g3_8 a3_8 b3_8 c4_4];
part7=[g3_8 c4_8 b3_8 d4_8 c4_4 c4_4];
para1=[part1 part2 part3 ];
para2=[part4 part5 ];
para3=[part6 part7 ];
legend=[para1 para2 para3];
sound(legend,fs)
将该法式在MATLAB中运行，我们可以听出音色一般，需要改良。
1.2  除噪音，加包络
下面经由过程加包络来消噪音。最简单的包络为指数衰减。最简单的指数衰减是对每个音乘以因子，在尝试中首先加的是的衰减，这种衰减方式利用的是不异速度的衰减，可是发现噪音并没有完全消弭，播放的音乐结果不是很好，感受音乐升沉性不强。于是采用分歧速度的衰减，按照乐音持续时候的长短来确定衰减的快慢，乐音持续时候越长，衰减的越慢，持续时候越短，衰减的越快。更科学的包络如下图所示，每个乐音都颠末冲激、衰减、持续、消逝四个阶段。
由上图可以看出这个包络是四段直线段组成的，是以只要确定了每段线段的端点，即可用端点数据写出直线方程，因为直线方程可以用通式写出（我用的是斜截式），是以这段包络可以用简单的轮回来完当作。例如认为包络线上的数据如下图所示：
据此在MATLAB中编写如下法式：
clear;clc;
fs=8000; %抽样频率
part1=[fre(55) fre(60) fre(60) fre(62) fre(64) fre(60) fre(64) fre(62) fre(59)];
part2=[fre(60) fre(60) fre(62) fre(64) fre(60) fre(59) fre(55)];
part3=[fre(60) fre(60) fre(62) fre(64) fre(65) fre(64) fre(62) fre(60) fre(59) fre(55) fre(57) fre(59) fre(60) fre(60)];
part4=[fre(57) fre(59) fre(57) fre(55) fre(57) fre(59) fre(60)];
part5=[fre(55) fre(57) fre(55) fre(53) fre(52) fre(53) fre(55)];
part6=[fre(57) fre(59) fre(57) fre(55) fre(57) fre(59) fre(60)];
part7=[fre(55) fre(60) fre(59) fre(62) fre(60) fre(60)];
para1=[part1 part2 part3 ];
para2=[part4 part5 ];
para3=[part6 part7 ];
f=[para1 para2 para3]; %各个乐音对应的频率
part1time=[0.5 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25];
part2time=[0.25 0.25 0.25 0.25 0.5 0.25 0.25];
part3time=[0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.5];
part4time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part5time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part6time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part7time=[0.25 0.25 0.25 0.25 0.5 0.5];
para1time=[part1time part2time part3time ];
para2time=[part4time part5time ];
para3time=[part6time part7time ];
time=fs*[para1time para2time para3time];%各个乐音的抽样点数
N=length(time); %这段音乐的总抽样点数
east=zeros(1,N); %用east标的目的量来储存抽样点
n=1;
for num=1:N %操纵轮回发生抽样数据，num暗示乐音编号
    t=1/fs:1/fs:time(num)/fs; %发生第num个乐音的抽样点
    baoluo=zeros(1,time(num)); %P为存储包络数据的标的目的量
for j=1:time(num)
if(j<0.2*time(num))
   y=7.5*j/time(num);
else
       if(j<0.333*time(num))
             y=-15/4*j/time(num)+9/4;
       else
if(j<0.666*time(num))
      y=1;
else
       y=-3*j/time(num)+3;
end
        end
end
baoluo(j)=y;
end
east(n:n+time(num)-1)=sin(2*pi*f(num)*t).*baoluo(1:time(num));
%给第num个乐音加上包络
    n=n+time(num);
end
sound(east,8000)       %播放音乐
plot(east)
此处操纵函数fre
function f = fre(p)
f=440*2^((p-69)/12);
运行获得的图像为：
下图是两个乐音交代处的局部放年夜图，可以看到知道前一个乐音衰减到零时，后一个乐音才起头从零增添，所以可以说消弭了噪音。
1.3  改编法式，实现1.2中乐曲的升八度
升高一个八度就是将每个乐音的频率都提高一倍，变为本来的2倍；降低一个八度即每个乐音的频率都减小一倍，变为本来的1/2。是以将存储乐音频率的标的目的量每个元素改变为2或1/2倍即可。
查表获得下列MATLAB法式
clc
clear
fs=44100;
t=0:1/fs:0.5;
c4_2=key(60, 2, fs);
d4_2=key(62, 2, fs);
e4_2=key(64, 2, fs);
f4_2=key(65, 2, fs);
g4_2=key(67, 2, fs);
a4_2=key(69, 2, fs);
b4_2=key(71, 2, fs);
c4_4=key(60, 4, fs);
d4_4=key(62, 4, fs);
e4_4=key(64, 4, fs);
f4_4=key(65, 4, fs);
g4_4=key(67, 4, fs);
a4_4=key(69, 4, fs);
b4_4=key(71, 4, fs);
c4_8=key(60, 8, fs);
d4_8=key(62, 8, fs);
e4_8=key(64, 8, fs);
f4_8=key(65, 8, fs);
g4_8=key(67, 8, fs);
a4_8=key(69, 8, fs);
b4_8=key(71, 8, fs);
c4_16=key(60, 16, fs);
d4_16=key(62,16, fs);
e4_16=key(64, 16, fs);
f4_16=key(65, 16, fs);
g4_16=key(67, 16, fs);
a4_16=key(69, 16, fs);
b4_16=key(71, 16, fs);
c5_2=key(72, 2, fs);
d5_2=key(74, 2, fs);
e5_2=key(76, 2, fs);
f5_2=key(77, 2, fs);
g5_2=key(79, 2, fs);
a5_2=key(81, 2, fs);
b5_2=key(83, 2, fs);
c5_4=key(72, 4, fs);
d5_4=key(74, 4, fs);
e5_4=key(76, 4, fs);
f5_4=key(77, 4, fs);
g5_4=key(79, 4, fs);
a5_4=key(81, 4, fs);
b5_4=key(83, 4, fs);
c5_8=key(72, 8, fs);
d5_8=key(74, 8, fs);
e5_8=key(76, 8, fs);
f5_8=key(77, 8, fs);
g5_8=key(79, 8, fs);
a5_8=key(81, 8, fs);
b5_8=key(83, 8, fs);
c5_16=key(72, 16, fs);
d5_16=key(74,16, fs);
e5_16=key(76, 16, fs);
f5_16=key(77, 16, fs);
g5_16=key(79, 16, fs);
a5_16=key(81, 16, fs);
b5_16=key(83, 16, fs);
part1=[g4_4 c5_8 c5_8 d5_8 e5_8 c5_8 e5_8 d5_8 b4_8];
part2=[c5_8 c5_8 d5_8 e5_8 c5_4 b4_8 g4_8 ];
part3=[c5_8 c5_8 d5_8 e5_8 f5_8 e5_8 d5_8 c5_8 b4_8 g4_8 a4_8 b4_8 c5_4 c5_4];
part4=[a4_8 b4_16 a4_8 g4_8 a4_8 b4_8 c5_4];
part5=[g4_8 a4_16 g4_8 f4_8 e4_8 f4_8 g4_4];
part6=[a4_8 b4_16 a4_8 g4_8 a4_8 b4_8 c5_4];
part7=[g4_8 c5_8 b4_8 d5_8 c5_4 c5_4];
para1=[part1 part2 part3 ];
para2=[part4 part5 ];
para3=[part6 part7 ];
legend=[para1 para2 para3];
sound(legend,fs)
1.4   在1.2的音乐中插手谐波
我们在1.2的音乐平分别加上二、三、四次谐波，其基波幅度为1，高次谐波幅度别离为0.2、0.3、0.1。应该将法式改为：
clear;clc;
fs=8000; %抽样频率
part1=[fre(55) fre(60) fre(60) fre(62) fre(64) fre(60) fre(64) fre(62) fre(59)];
part2=[fre(60) fre(60) fre(62) fre(64) fre(60) fre(59) fre(55)];
part3=[fre(60) fre(60) fre(62) fre(64) fre(65) fre(64) fre(62) fre(60) fre(59) fre(55) fre(57) fre(59) fre(60) fre(60)];
part4=[fre(57) fre(59) fre(57) fre(55) fre(57) fre(59) fre(60)];
part5=[fre(55) fre(57) fre(55) fre(53) fre(52) fre(53) fre(55)];
part6=[fre(57) fre(59) fre(57) fre(55) fre(57) fre(59) fre(60)];
part7=[fre(55) fre(60) fre(59) fre(62) fre(60) fre(60)];
para1=[part1 part2 part3 ];
para2=[part4 part5 ];
para3=[part6 part7 ];
f=[para1 para2 para3]; %各个乐音对应的频率
part1time=[0.5 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25];
part2time=[0.25 0.25 0.25 0.25 0.5 0.25 0.25];
part3time=[0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.5];
part4time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part5time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part6time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part7time=[0.25 0.25 0.25 0.25 0.5 0.5];
para1time=[part1time part2time part3time ];
para2time=[part4time part5time ];
para3time=[part6time part7time ];
time=fs*[para1time para2time para3time];%各个乐音的抽样点数
N=length(time); %这段音乐的总抽样点数
east=zeros(1,N); %用east标的目的量来储存抽样点
n=1;
for num=1:N %操纵轮回发生抽样数据，num暗示乐音编号
    t=1/fs:1/fs:time(num)/fs; %发生第num个乐音的抽样点
    baoluo=zeros(1,time(num)); %P为存储包络数据的标的目的量
for j=1:time(num)
if(j<0.2*time(num))
   y=7.5*j/time(num);
else
       if(j<0.333*time(num))
             y=-15/4*j/time(num)+9/4;
       else
if(j<0.666*time(num))
      y=1;
else
       y=-3*j/time(num)+3;
end
        end
end
baoluo(j)=y;
end
         h=[1 0.2 0.3 0.1];
         xiebo=zeros(1,length(t));
for i=1:length(n)
      xiebo=xiebo+h(i)*sin(2*i*pi*f(num)*t);
end
east(n:n+time(num)-1)=xiebo.*baoluo(1:time(num));
n=n+time(num);
end
sound(east,8000) %播放音乐
plot(east) %图像
插手谐波后音色获得较着好转
运行获得的图像为：
2
2.用傅里叶变换阐发音乐
2.1 阐发can.wav的调子和节奏
我们对can进行傅里叶变换阐发其基波协调波，获得can的幅值谱，频谱图上的第一个凸起的波峰对应的频率即为can的基频，可编写了如下法式：
clear;clc;
[y,Fs]= audioread('can.wav');
fs=8000;
NFFT = 2^nextpow2(length(y));
Y = fft(y,NFFT)/length(y);
g = fs/2*linspace(0,1,NFFT/2+1);
plot(g,2*abs(Y(1:NFFT/2+1)))
运行后获得的成果为
经由过程增添can的周期性显示出离散化水平高的幅值谱，即让can在时域反复多次后在进行傅里叶变换。
操纵repmat函数将can在时域内反复。法式可点窜为:
clear;clc;
WAV= audioread('can.wav');
fs=8000;
wave2proc =repmat(WAV,20,1); %将 can反复20次
NFFT = 2^nextpow2(length(WAV));
Y = fft(WAV,NFFT)/length(WAV);
g = fs/2*linspace(0,1,NFFT/2+1);
plot(g,2*abs(Y(1:NFFT/2+1)))
由图读出can的基频为329.1Hz，幅值为2.451，高次谐波幅值别离为：
2.2 按照快速傅里叶变换合当作音乐
将法式中的波形幅度矩阵
m=[1 0.3 0.2 0.1]
改为
m=[2.0912 4.0597 4.7156 7.5215 5.6484 4.9845 3.4894 2.4568];
即可
法式如下
clear;clc;
fs=8000; %抽样频率
part1=[fre(67) fre(72) fre(72) fre(74) fre(76) fre(72) fre(76) fre(74) fre(71)];
part2=[fre(72) fre(72) fre(74) fre(76) fre(72) fre(71) fre(67)];
part3=[fre(72) fre(72) fre(74) fre(76) fre(77) fre(76) fre(74) fre(72) fre(71) fre(67) fre(69) fre(71) fre(72) fre(72)];
part4=[fre(69) fre(71) fre(69) fre(67) fre(69) fre(71) fre(72)];
part5=[fre(67) fre(69) fre(67) fre(65) fre(64) fre(65) fre(67)];
part6=[fre(69) fre(71) fre(69) fre(67) fre(69) fre(71) fre(72)];
part7=[fre(67) fre(72) fre(71) fre(74) fre(72) fre(72)];
para1=[part1 part2 part3 ];
para2=[part4 part5 ];
para3=[part6 part7 ];
f=[para1 para2 para3]; %各个乐音对应的频率
part1time=[0.5 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25];
part2time=[0.25 0.25 0.25 0.25 0.5 0.25 0.25];
part3time=[0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.5];
part4time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part5time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part6time=[0.25 0.125 0.25 0.25 0.25 0.25 0.5];
part7time=[0.25 0.25 0.25 0.25 0.5 0.5];
para1time=[part1time part2time part3time ];
para2time=[part4time part5time ];
para3time=[part6time part7time ];
time=fs*[para1time para2time para3time];%各个乐音的抽样点数
N=length(time); %这段音乐的总抽样点数
east=zeros(1,N); %用east标的目的量来储存抽样点
n=1;
for num=1:N %操纵轮回发生抽样数据，num暗示乐音编号
    t=1/fs:1/fs:time(num)/fs; %发生第num个乐音的抽样点
    baoluo=zeros(1,time(num)); %P为存储包络数据的标的目的量
for j=1:time(num)
if(j<0.2*time(num))
   y=7.5*j/time(num);
else
       if(j<0.333*time(num))
             y=-15/4*j/time(num)+9/4;
       else
if(j<0.666*time(num))
      y=1;
else
       y=-3*j/time(num)+3;
end
        end
end
baoluo(j)=y;
end
         h=[2.0912 4.0597 4.7156 7.5215 5.6484 4.9845 3.4894 2.4568];
         xiebo=zeros(1,length(t));
for i=1:length(n)
      xiebo=xiebo+h(i)*sin(2*i*pi*f(num)*t);
end
east(n:n+time(num)-1)=xiebo.*baoluo(1:time(num));
n=n+time(num);
end
sound(east,8000) %播放音乐
plot(east) %图像
3
第三部门音乐旌旗灯号的处置
3.1 插手延时和混响
选择一段语音旌旗灯号作为阐发对象，并对其进行频谱阐发，在时域顶用数字旌旗灯号处置方式给旌旗灯号插手3种混响，再阐发其频谱，并与原始旌旗灯号进行比力。
设计思绪
1、操纵Windows下的灌音机或其他软件，录制一段语音旌旗灯号，时候节制在3s摆布，并对录制的旌旗灯号进行采样
2、语音旌旗灯号的频谱阐发，画出采样后的时域波形和频谱图
3、将旌旗灯号插手延时和混响，再阐发当时域波形和频谱图，并与原始旌旗灯号频谱进行比力
所以按照设计思绪进行尝试
1、读取3s的语音旌旗灯号并画出时域波形和频谱图
x1=audioread('good.wav');
[x,fs]=audioread('good.wav');
x=x(:,1);   %取单声道
sound(x,fs);
X=fft(x,640000); %gaidong
magX=abs(X);
angX=angle(X);
figure(1);
subplot(2,1,1);plot(x);title('yuanshi boxing');
subplot(2,1,2);plot(X);title('yuanshi pinpu');
获得的波形和频谱图为
2、对语音旌旗灯号进行采样并画出采样后旌旗灯号的时域波和频谱图
[x,fs]=audioread('good.wav');
x=x(:,1);
sound(5*x,fs);
n1=0:2000;
N=size(x,1);
Y=fft(x,320000);
figure(2);
subplot(2,1,1);plot(x);title('caiyang boxing');
subplot(2,1,2);plot(n1(1:1000),Y(1:1000));title('caiyang pinpu');
获得的采样波形和频谱为
3、对采样后的旌旗灯号延时，并画出延时后的时域波形和频谱图
z1=[zeros(10000,1);x];           %对旌旗灯号进行延时
z2=[zeros(20000,1);x];
z3=[zeros(30000,1);x];
Z1=fft(z1,160000);
Z2=fft(z2,160000);
Z3=fft(z3,160000);
figure(3);
subplot(3,1,1);plot(z1); title('延时后的时域图1');        %画出延时后的旌旗灯号时域图
subplot(3,1,2);plot(z2); title('延时后的时域图2');
subplot(3,1,3);plot(z3); title('延时后的时域图3');
figure(4)
subplot(3,1,1);plot(n1(1:1000),Z1(1:1000));title('延时后的频谱图1');     %延时后的旌旗灯号频谱图
subplot(3,1,2);plot(n1(1:1000),Z2(1:1000));title('延时后的频谱图2');
subplot(3,1,3);plot(n1(1:1000),Z3(1:1000));title('延时后的频谱图3');
获得延时的时域图和频谱图为
4、对旌旗灯号进行混响，并画出混响后的时域波形和频谱图
z1=[zeros(10000,1);x];           %对旌旗灯号进行延时
z2=[zeros(20000,1);x];
z3=[zeros(30000,1);x];
x1=[x;zeros(10000,1)];                 %使语音旌旗灯号与延时后旌旗灯号划一长度
x2=[x;zeros(20000,1)];
x3=[x;zeros(30000,1)];
y1=plus(x1,z1);                             %旌旗灯号的混响
y2=plus(x2,z2);
y3=plus(x3,z3);
sound(y1,fs);
sound(y2,fs);
sound(y3,fs);
figure(5);
subplot(3,1,1);plot(y1); title('混响的时域图1');     %混响时域图
subplot(3,1,2);plot(y2); title('混响的时域图2');
subplot(3,1,3);plot(y3); title('混响的时域图3');
Y1=fft(y1,160000);                                %对混响旌旗灯号FFT变换
Y2=fft(y2,160000);
Y3=fft(y3,160000);
figure(6);
subplot(3,1,1);plot(n1(1:1000),Y1(1:1000)); title('混响的频谱图1');    %混响频谱图
subplot(3,1,2);plot(n1(1:1000),Y2(1:1000)); title('混响的频谱图2');
subplot(3,1,3);plot(n1(1:1000),Y3(1:1000)); title('混响的频谱图3');
获得的混响后的时域波形和频谱图为
3.2 平衡处置
设计思绪：设计三个混响器作为平衡处置的东西
（1）无限反响混响器
a=0.05;                               %a取小于等于1
Bz=[0,0,0,0,0,0,0,0,0,0,1];            %分子的系数
Az=[1,0,0,0,0,0,0,0,0,0,-a];           %分母的系数
yy1=filter(Bz,Az,z1);                  %滤波器进行滤波
YY1=fft(yy1,320000);    %经无限反响滤波器后的旌旗灯号做32000点的FFT变换
（2）多重反响混响器
a=0.05;                            %a取小于等于1
N=5
Bz1=[1,0,0,0,0,0,0,0,0,0,-0.5^N]      %分子的系数
Az1=[1,0,0,0,0,0,0,0,0,0,-0.5];        %分母的系数
yy2=filter(Bz1,Az1,z1);                 %滤波器进行滤波
YY2=fft(yy2,320000);    %经多重反响滤波器后的旌旗灯号做32000点的FFT变换
（3）全通布局的混响器
a=0.05;                               %a取小于等于1
Bz1=[a,0,0,0,0,0,0,0,0,0,1];          %分子的系数
Az1=[1,0,0,0,0,0,0,0,0,0,a];          %分母的系数
yy3=filter(Bz1,Az1,z1);               %滤波器进行滤波
YY3=fft(yy3,320000);    %经全通布局的混响器后的旌旗灯号做32000点的FFT变
最后输作声音
sound(yy1,fs);
sound(yy2,fs);
sound(yy3,fs);
（4）画出经混响器处置后旌旗灯号的时域波形和频谱图
figure(8);
subplot(2,1,1);
plot(yy1);
title('无限个反响滤波器时域图');       %无限反响滤波器时域波形
subplot(2,1,2);
plot(n1(1:1000),YY1(1:1000));
title('无限个反响滤波器频谱图 ');    %无限反响滤波器频谱图
figure(9)
subplot(2,1,1);
plot(yy2);
title('多重反响滤波器的时域图')       %多重反响滤波器的混响器时域波形
subplot(2,1,2);
plot(n1(1:1000),YY2(1:1000));
title('多重反响滤波器的频谱图')    %多重反响滤波器的频谱图
figure(10)
subplot(2,1,1);
plot(yy3);
title('全通布局滤波器的时域图')     %全通布局的混响器时域波形
subplot(2,1,2);
plot(n1(1:1000),YY3(1:1000));
title('全通布局滤波器的频谱图')     %全通布局的混响器频谱图
获得的波形为