MATLAB语音信号处理

数字信号处理课设，我们使用MATLAB对语音信号进行了一系列处理，并将其所有功能集中于下图界面中：

这个界面涉及功能众多，其中包括语音信号的观察分析、音色变换、AM调制解调、减抽样、加噪去噪、相频分析和幅频滤波等，最重要的是对MATLAB中函数的掌握，通过不同函数的组合实现你想要实现的功能。

本篇不会给出整个界面的程序，下面会分块给出每个功能的程序，整个界面只需GUI设计界面文件、定义结构体并把对应键程序打进去即可。

1、语音信号的采集

1.1题目要求

使用windows下的录音机录制一段语音信号、音乐信号或者采用其他软件截取一段音乐信号（要求：时间不超过5s，文件格式为WAV。）

① 请每位同学都参与录音，内容自定。

② 使用wavread语句读取语音/音乐信号获取抽样率；（注意：读取的信号是双声道信号，即为双列向量，需要分列处理）；

③ 输出时域语音/音乐信号的波形。

④ 实现对录音信号的声音大小的调节。

⑤ 实现对两种语音/音乐信号的混音音效。

⑥ 实现音乐信号的回音音效。

1.2设计内容及方案

① 读取音频信号：我是通过wavread函数读取.wav文件的方式来获得，当然首先要自己创建一个.wav音频，我是通过电脑录音生成.mp3然后格式工厂转成.wav的，需保存到同一文件夹下。

② 分声道处理：一般音乐和语音信号都是双声道信号，时域和频谱图会有两个颜色，所以要取单列来分析，通过x1=x(:,1)语句来实现。

③ 画时域波形图：用plot函数来画图，注意横坐标为时间t。

④ 音量大小调节：通过将音频直接乘一个系数来实现调音量。

⑤ 混音和回声：混音即将两个音频相加，要相加就得保证矩阵一样，所以要通过截取并补零矩阵来实现；回声是把三个信号叠加，这三个信号在不同位置补零音量也逐渐变小，就可以实现回声。

⑥ 播放声音：本题我使用wavplay来播放声音，会有警告，后面的题我用sound比较好。

1.3程序源码及注释

clear
[x,fs] = wavread('beautiful.wav');%音乐信号
[y,fs1]= wavread('1.wav');%女生声音
[z,fs2]= wavread('2.wav');%男生声音
%输出频率
fs
fs1
fs2
%音乐语音信号分声道处理
x1=x(:,1);
y1=y(:,1);
z1=z(:,1);
%画音乐信号时域图
n1=length(x1);%length取数列长度即元素个数
figure(1)
t1=(0:(n1-1))/fs;
plot(t1,x1);
axis([0,5,-1,1]);
xlabel('时间t');
ylabel('幅度');
title('音乐信号时域波形');
%画语音信号时域图
n2=length(y1);
figure(2)
subplot(2,1,1);
t2=(0:(n2-1))/fs1;
plot(t2,y1);%女生
axis([0,4,-0.5,0.5]);
xlabel('时间t');
ylabel('幅度');
title('女生语音信号时域波形');
n3=length(z1);
subplot(2,1,2);
t3=(0:(n3-1))/fs2;
plot(t3,z1);%男生
axis([0,4,-0.5,0.5]);
xlabel('时间t');
ylabel('幅度');
title('男生语音信号时域波形');
%对语音信号声音大小调节
wavplay(y,fs1); %播放原语音
y11=10*y;
wavplay(y11,fs1); %加大音量播放
y22=0.5*y;
wavplay(y22,fs1); %减小音量播放
%两种语音信号的混音
[m,n]=size(y1);%size取矩阵的行列数
[m0,n0]=size(z1);
a=zeros(abs(m-m0),n);%两矩阵行数差为零矩阵行数
if length(y1)<length(z1)
     y2=[y1;a];
     y3=y2+z1;%两个矩阵行数一样才能相加，所以要补零
else 
     y2=[z1;a];
     y3=y2+y1;%y1和z1中长的那个不变，短的那个补零
end;
wavplay(y3,fs1) ;%播放混音语音
  %画混音波形
figure(3)
subplot(2,1,1);
t4=(0:(max(n2,n3)-1))/fs1;
plot(t4,y3);
axis([0,4.5,-0.5,0.5]);
xlabel('时间');
ylabel('幅度');
title('两语音信号叠加后时域波形');
%音乐信号的回音
x11=x1(1:200000);%截取部分
x11=x11';%因为输出为一列所以要转置成一行
int0=zeros(1,20000);%1行2000列的零矩阵
temp1=[x11,int0,int0];
temp2=[int0,0.6*x11,int0];
temp3=[int0,int0,0.3*x11];%通过补零实现延时，同时声音一个比一个小
hui=temp1+temp2+temp3;%三重声音相加实现回声
N=length(hui);
wavplay(hui,fs1);%播放回音音乐
  %画回声波形
subplot(2,1,2);
t1=(0:(N-1))/fs;
plot(t1,hui);
axis([0,4.5,-1,1]);
xlabel('时间');
ylabel('幅度');
title('回声时域波形');

1.4运行结果

仿真结果分析：我听到了原声和音量放大减小的声音，也听到了混音和回声的效果，变化明显；本题我画了音乐和两个语音信号的时域波形以及混音回声的时域波形，音乐信号幅度比语音信号高且连贯性高，混音之后幅度叠加，回声之后幅度也增大，波形有很明显的变化。

2、语音/音乐信号的频谱和音谱的观察

2.1题目要求

① 输出语音/音乐信号的波形和频谱，观察现象；

② 使用sound语句播放语音/音乐信号，注意不同抽样率下音调变化，解释现象。

2.2设计内容及方案

本题读取音频信号、画时域波形和播放原理和上题一样，涉及的新内容有：

① 画频谱图：我将横坐标设为频率f，纵坐标需要用fft函数求傅里叶变换然后利用abs函数求幅值画幅度谱，再用plot画出频谱图。

② 调节频率：我将频率fs乘一个系数放大缩小并播放，感受频率对语音的影响。

2.3程序源码及注释

clear
[x,fs] = wavread('beautiful.wav');%音乐信号
[y,fs1]= wavread('1.wav');%女生声音
[z,fs2]= wavread('2.wav');%男生声音
%输出频率
fs
fs1
fs2
%音乐语音信号分声道处理
x1=x(:,1);
y1=y(:,1);
z1=z(:,1);
%画音乐信号时域图
n1=length(x1);%length取数列长度即元素个数
figure(1)
subplot(2,1,1);
t1=(0:(n1-1))/fs;
plot(t1,x1);
axis([0,5,-1,1]);
xlabel('时间t');
ylabel('幅度');
title('音乐信号时域波形');
%画音乐信号频域图
X1=fft(x1,n1);
subplot(2,1,2);
f1=0:fs/n1:fs*(n1-1)/n1;
plot(f1,abs(X1));%也可以使用fftshift函数使fft后的结果以fs/2为中心左右互换
axis([0,44100,0,6000]);
xlabel('频率f');
ylabel('幅度');
title('音乐信号频谱');
%画女生语音信号时域图
n2=length(y1);
figure(2)
subplot(2,2,1);
t2=(0:(n2-1))/fs1;
plot(t2,y1);
axis([0,4,-0.5,0.5]);
xlabel('时间t');
ylabel('幅度');
title('女生语音信号时域波形');
%画女生语音信号频域图
Y=fft(y1,n2);
subplot(2,2,2);
f2=0:fs1/n2:fs1*(n2-1)/n2;
plot(f2,abs(Y));
axis([0,48000,0,700]);
xlabel('频率f');
ylabel('幅度');
title('女生语音信号频谱');
%画男生语音信号时域图
n3=length(z1);
subplot(2,2,3);
t3=(0:(n3-1))/fs2;
plot(t3,z1);
axis([0,4,-0.5,0.5]);
xlabel('时间t');
ylabel('幅度');
title('男生语音信号时域波形');
%画男生语音信号频域图
Z=fft(z1,n3);
subplot(2,2,4);
f3=0:fs2/n3:fs2*(n3-1)/n3;
plot(f3,abs(Z));
axis([0,48000,0,700]);
xlabel('频率f');
ylabel('幅度');
title('男生语音信号频谱');
%不同抽样率语音信号
sound(y1,fs);%播放原音乐
pause(5);
sound(y1,2*fs);%播放高频率音乐
pause(2);
sound(y1,0.5*fs);%播放低频率音乐

2.4运行结果

仿真结果分析：本题中画了频谱图，可以更直观的看到信号的特征，包括截至频率和音频范围，我的音乐和语音信号主要集中在低频段，而且音乐信号的幅度要比语音信号的幅度高。通过放大缩小频率的声音变化，即频率大时语调变高且语速加快，让我感受到了频率对信号的作用。

3、语音/音乐信号的抽取（减抽样）

3.1题目要求

① 观察语音/音乐信号频率的上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；

② 输出减抽样语音/音乐信号的波形和频谱，观察现象，给出理论解释；

③ 播放减抽样语音/音乐信号，注意抽样率改变，比较不同抽取间隔下的声音，解释现象。

3.2设计内容及方案

减抽样：本题研究的唯一内容就是对信号以不同间隔抽样，我抽取了200000长度的信号，便于以一定间隔抽样。首先我以小间隔2抽样，信号声音基本和原声差不多，没有发生混叠；而后我又用大间隔20抽样，信号声音有了很明显的变化，即发生了混叠。

3.3程序源码及注释

clear
[x,fs]=wavread('beautiful.wav');%音乐信号
x1=x(:,1);%取单列
x2=x1(1:200000);%截取长度
N1=length(x2)%求信号长度
%画原信号时域波形和频谱
t1=(0:(N1-1))/fs;
figure(1);
subplot(2,1,1);
plot(t1,x2);
title('原信号时域波形');
xlabel('时间t');
ylabel('幅度');
f1=0:fs/N1:fs*(N1-1)/N1;
Fx1=fft(x2,N1); %求傅里叶变换
subplot(2,1,2);
plot(f1,abs(Fx1));
title('原信号的频谱');%画原信号频谱
xlabel('频率f');
ylabel('幅度');
%非混叠间隔减抽样
d1=2;j=0;
for i=1:d1:200000
    j=j+1;
    x22(j)=x2(i);%对信号以d1为间隔做减抽样
end
N2=length(x22);%对信号以d1为间隔做减抽样的长度
  %画d1间隔抽样图
t2=(0:(N2-1))/fs;
figure(2);
subplot(2,1,1);
plot(t2,x22);
title('以d1为间隔减抽样时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,1,2);
f2=0:(fs/d1)/N2:(fs/d1)*(N2-1)/N2;
Fx2=fft(x22,N2);%求傅里叶变换画频谱
plot(f2,abs(Fx2));
title('以d1为间隔减抽样后的频谱(非混叠)');
xlabel('频率f');
ylabel('幅度');
%混叠间隔减抽样
d2=20;j=0;
for i=1:d2:200000
    j=j+1;
    x3(j)=x2(i);%对信号以d2为间隔做减抽样
end
N3=length(x3) ;%对信号以d2为间隔做减抽样长度
  %画d2间隔抽样图
t3=(0:(N3-1))/fs;
figure(3)
subplot(2,1,1);
plot(t3,x3);
title('以d2为间隔减抽样后的时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,1,2);
f3=0:(fs/d2)/N3:(fs/d2)*(N3-1)/N3;
Fx3=fft(x3,N3);%求傅里叶变换画频谱
plot(f3,abs(Fx3));
title('以d2为间隔减抽样后的频谱(混叠)');
xlabel('频率f');
ylabel('幅度');
%播放减抽样信号
sound(x2,fs);%原信号播放
pause(5);
sound(x22,fs/d1);%以d1为间隔减抽样信号播放(非混叠)
pause(8);
sound(x3,fs*(1/d2));%以d2为间隔减抽样信号播放(混叠)

3.4运行结果

仿真结果分析：从信号的时域波形和频谱图可以看出抽样后时间范围和频率范围都有所减小，小间隔抽样频谱波形变化比较小，没有发生混叠，大间隔抽样频谱波形有了很大的变化，发生了混叠。抽取间隔越小，声音越清晰，时间间隔越大，声音越不清晰，混叠现象越明显。未混叠时，声音尖锐，混叠时，声音轻，只有淡淡的音调，基本没有起伏，不清晰。

原因：采样频率fs必须大于等于2fc，采样频率小于2fc时，发生混叠，若采样频率越小，则混叠越明显。我的图中fc不到5000Hz,间隔2采样频率fs为22000Hz大于两倍的fc所以不混叠；间隔20采样频率fs为2200Hz小于两倍的fc所以发生混叠。

4 语音/音乐信号的AM调制

4.1题目要求

① 观察语音/音乐信号的频率上限，选择适当调制频率对信号进行调制（给出高、低两种调制频率）；

② 输出调制信号的波形和频谱，观察现象，给出理论解释；

③ 播放调制语音/音乐信号，注意不同调制频率下的声音，解释现象。

4.2设计内容及方案

在这道题中我统一使用了频率归一化，便于从原信号中读取截止频率和设置载波频率。

① 取低频载波对信号进行AM调制：这里取了0.1pi为低频调制载波频率，与原信号相乘实现AM调制，这里用点乘转置矩阵实现。

② 取高频载波对信号进行AM调制：这里取了0.7pi为低频调制载波频率，与原信号相乘实现AM调制，这里用点乘转置矩阵实现。

③ 播放调制后信号：分别播放低频和高频调制时的音乐，用sound函数播放。

4.3程序源码及注释

clear
[x,fs]=wavread('beautiful.wav'); %读取音乐信号
x1=x(:,1);%取单列 
N=length(x1);%求信号长度 
n=0:N-1;%所有元素
t=(0:(N-1))/fs;%时间
f=0:fs/N:fs*(N-1)/N;%频率
w=2*f/fs;%归一化
%画原信号时域波形和频谱
figure(1);
subplot(2,1,1);
plot(t,x1);
title('原信号时域波形');
xlabel('时间t');
ylabel('幅度');
Fx1=fft(x1,N); %求傅里叶变换
subplot(2,1,2);
plot(w,abs(Fx1));
title('原信号的频谱');%画原信号频谱
xlabel('w');
ylabel('幅度');
%低频调制
x2=cos(n*0.1*pi);%低频时余弦信号
x22=x1.*x2';%低频调制信号
X2=fft(x2); 
X22=fft(x22);
%高频调制
x3=cos(n*0.7*pi);%高频时余弦信号
x33=x1.*x3';%高频调制信号
X3=fft(x3); 
X33=fft(x33); 
%画出低频时载波和调制信号波形
figure(2)
subplot(2,2,1);
plot(t,x2);
axis([0,0.001,-1,1]);
title('低频时余弦信号时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,2,2);
plot(w,abs(X2));
title('低频时余弦信号频谱');
xlabel('w');
ylabel('幅度');

subplot(2,2,3);
plot(t,x22);
title('低频调制信号时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,2,4);
plot(w,abs(X22));
title('低频调制信号频谱');
xlabel('w');
ylabel('幅度');
%画出高频时载波和调制信号波形
figure(3)
subplot(2,2,1);
plot(t,x3);
axis([0,0.001,-1,1]);
title('高频时余弦信号时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,2,2);
plot(w,abs(X3));
title('高频时余弦信号频谱');
xlabel('w');
ylabel('幅度');

subplot(2,2,3);
plot(t,x33);
title('高频调制信号时域波形');
xlabel('时间t');
ylabel('幅度');
subplot(2,2,4);
plot(w,abs(X33));
title('高频调制信号频谱');
xlabel('w');
ylabel('幅度');
%播放调制信号
sound(x1,fs);%播放原音乐
pause(5);
sound(x22,fs);%低频调制播放
pause(5);
sound(x33,fs);%高频调制播放

4.4运行结果

仿真结果分析：从图中可以看出原音乐信号的截止频率为0.2pi，这里设置了0.1pi和0.7pi两种频率对信号进行AM调制，原信号的调制相当于频谱搬移， 左移一个右移一个，调制的目的是便于信号在信道中传输。当调制频率较高时，声音响度低，几乎只能听见兹兹的声音， 信号几乎完全失真，当调制频率较低时，声音很尖锐，响度较大，能听出调子，但也有兹兹的声音。

5、AM调制语音/音乐信号的同步解调

5.1题目要求

① 设计巴特沃斯滤波器完成同步解调，观察滤波器频率响应曲线；

② 窗函数法设计FIR滤波器完成同步解调，观察滤波器频率响应曲线（要求：分别使用矩形窗和布莱克曼窗，进行比较）；

③ 输出解调信号的波形和频谱，观察现象，给出理论解释；

④ 播放解调语音/音乐信号，比较不同滤波器下的声音，解释现象。

5.2设计内容及方案

① 对调制后的信号进行解调：将调制后的信号与调制时相同的载波相乘实现解调，这里用点乘转置矩阵实现。

② 用巴特沃斯滤波器对解调信号进行滤波：首先求巴特沃斯滤波器的频率响应，其中用到了buttord求满足性能指标的滤波器阶数N和3dB截止频率wc、用butter计算模拟滤波器的传输函数Ha(s)、用freqz求频响。然后用filter实现滤波。

③ 用矩形窗FIR滤波器对解调信号进行滤波：首先求矩形窗FIR滤波器的频率响应，其中先求理想低通单位脉冲响应hd，然后加矩形窗截断求模拟脉冲响应，再利用freqz求频率响应。然后利用卷积conv实现滤波。

④ 用布莱克曼窗FIR滤波器对解调信号进行滤波：首先求布莱克曼窗FIR滤波器的频率响应，其中先求理想低通单位脉冲响应hd，然后加布莱克曼窗截断求模拟脉冲响应，再利用freqz求频率响应。然后利用卷积conv实现滤波。

⑤ 播放调制、解调和滤波后的声音：通过声音变化感受调制、解调和滤波。

5.3程序源码及注释

clear
[y,fs]=audioread('beautiful.wav');%读取音乐信号
y1=y(:,1);%取单列 
N1=length(y1);%求信号长度 
n=0:N1-1;%所有元素
t=n/fs;%时间
f=0:fs/N1:fs*(N1-1)/N1;%频率
w=2*f/fs;%归一化
y2=cos(n*pi*0.12);%0.12pi时余弦信号  
Y=y1.*y2';%音乐信号AM调制
Y2=Y.*y2';%解调相乘器
Y22=fft(Y2);%傅里叶变换
%画解调信号时域波形和频谱
figure(1);
subplot(2,1,1);
plot(t,Y2);
xlabel('时间t');
ylabel('幅度');
title('解调信号时域波形');
subplot(2,1,2);
plot(w,abs(Y22));
xlabel('w');
ylabel('幅度');
title('解调信号频谱');
%求巴特沃斯滤波器频率响应
wp=0.12;ws=0.3;rp=1;rs=50; %设计巴特沃斯IIR滤波器参数
[N,wc]=buttord(wp,ws,rp,rs); %求满足性能指标的滤波器阶数N和3dB截止频率wc
[b,a]=butter(N,wc); %计算模拟滤波器的传输函数Ha(s)
[Hd,w]=freqz(b,a);%求频响
figure(2)
subplot(3,1,1);
plot(w/pi,abs(Hd));
axis([0,1,0,1.5]);
xlabel('w/π');
ylabel('幅度');
title('巴特沃斯频率响应曲线');
%求巴特沃斯滤波器滤波信号
Y3=filter(b,a,Y2);%滤波音乐信号
N3=length(Y3)
t1=(0:(N3-1))/fs;%时间
f1=0:fs/N3:fs*(N3-1)/N3;%频率
w1=2*f1/fs;%归一化
Y33=fft(Y3);%傅里叶变换
subplot(3,1,2);
plot(t1,Y3);
xlabel('时间t');
ylabel('幅度');
title('巴特沃斯滤波信号时域波形');
subplot(3,1,3);
plot(w1,abs(Y33));
xlabel('w');
ylabel('幅度');
title('巴特沃斯滤波信号频谱');
%理想低通单位脉冲响应hd计算(窗函数法要用)
N=11;n1=[0:1:(N-1)];wc1=0.2*pi;  
alpha=(N-1)/2;
m=n1-alpha+eps;%加一个小数避免0作除数
hd=sin(wc1*m)./(pi*m);
%求矩形窗FIR滤波器频率响应
w_han1=(boxcar(N))';%矩形窗
h1=hd.*w_han1;%加窗截断得脉冲响应
    %以下为频率响应计算
[H,w]=freqz(h1,1,1000,'whole');
mag=abs(H);%mag为幅度
db=20*log10((mag+eps)/max(mag));
pha=angle(H);%pha为相位
grd=grpdelay(h1,1,w);%grd为群延迟
    %以下为画图
figure(3);
subplot(3,2,1);
plot(w/pi,db);
axis([0,2,-200,30]);
xlabel('w/π');
ylabel('幅度');
title('矩形窗频响图');
%求矩形窗FIR滤波器滤波信号
Y4=conv(Y2,h1);%卷积滤波
N4=length(Y4);
t2=(0:(N4-1))/fs;%时间
f2=0:fs/N4:fs*(N4-1) /N4;%频率
w2=2*f2/fs;%归一化
Y44=fft(Y4,N4);%傅里叶变换
subplot(3,2,3);
plot(t2,Y4);
xlabel('时间t');
ylabel('幅度');
title('矩形窗滤波信号时域波形');
subplot(3,2,5);
plot(w2,abs(Y44));
xlabel('w');
ylabel('幅度');
title('矩形窗滤波信号频谱');
%求布莱克曼窗FIR滤波器频率响应
w_han2=(blackman(N))';%布莱克窗
h2=hd.*w_han2;%加窗截断得脉冲响应
     %以下为频率响应计算
[H,w]=freqz(h2,1,1000,'whole');
mag=abs(H);%mag为幅度
db=20*log10(mag+eps)/max(mag);
pha=angle(H);%pha为相位
grd=grpdelay(h2,1,w);%grd为群延迟
     %以下为画图
subplot(3,2,2);
plot(w/pi,db);
axis([0,2,-200,30]);
xlabel('w/π');
ylabel('幅度');
title('布莱克曼窗频响图');
%求布莱克曼窗FIR滤波器滤波信号
Y5=conv(Y2,h2);%卷积滤波
N5=length(Y5);
t3=(0:(N4-1))/fs;%时间
f3=0:fs/N4:fs*(N4-1)/N4;%频率
w3=2*f3/fs;%归一化
Y55=fft(Y5,N5);%傅里叶变换
subplot(3,2,4);
plot(t3,Y5);
xlabel('时间t');
ylabel('幅度');
title('布莱克曼滤波信号时域波形');
subplot(3,2,6);
plot(w3,abs(Y55));
xlabel('w');
ylabel('幅度');
title('布莱克曼滤波信号频谱图');
%播放解调后声音
sound(y1,fs);%原声音
pause(5);
sound(Y,fs);%AM调制后的声音 
pause(5);
sound(Y2,fs);%解调后的声音 
pause(5);
sound(Y3,fs);%巴特沃斯滤波器滤波后的声音 
pause(5);
sound(Y4,fs);%矩形窗滤波器滤波后的声音 
pause(5);
sound(Y5,fs);%布莱克曼窗滤波后的声音

5.4运行结果

 仿真结果分析：首先得到了解调后的时域波形和频谱图，可以看出解调后的信号并没有完全恢复原信号，会夹杂一点调制过程中的载波，通过滤波后信号频谱有了很大改善。播放声音发现：巴特沃斯滤波后声音清晰，基本和原来的音乐差不多，但是音乐稍微低沉。巴特沃斯滤波器的特点是通频带的频率响应曲线平滑。矩形窗滤波声音较为沉闷，也伴有杂音。原因是矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低 。布莱克曼窗滤波声音与原信号较为接近，声音有些尖锐，也有轻微杂音。原因是布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。

6、语音/音乐信号的滤波去噪

6.1题目要求

① 原始信号叠加幅度为0.05，频率为3kHz,5kHz,8kHz的三余弦混合噪声，观察噪声频谱以及加噪后语音/音乐信号的音频和频谱，并播放音乐，感受噪声对语音/音乐信号的影响；

② 给原始语音/音乐信号叠加幅度为0.5的随机白噪声（可用rand语句产生），观察噪声频谱以及加噪后语音/音乐信号的音频和频谱，并播放音乐，感受噪声对语音/音乐信号的影响；

③ 根据步骤①、②观察到的频谱，选择合适指标设计滤波器进行滤波去噪，关注去噪后信号音谱和频谱，并播放音乐，解释现象。

6.2设计内容及方案

① 给信号加三余弦混合噪声：首先求得三余弦混合噪声，幅值为0.05，分别给频率3kHz,5kHz,8kHz，叠加得到噪声，然后将噪声和原信号叠加。

② 给信号加随机白噪声：首先求得随机白噪声，幅值为0.5，用randn语句产生噪声，然后将噪声和原信号叠加。

③ 滤掉噪声：我使用了巴特沃斯滤波器来滤噪，其中用到buttord求满足性能指标的滤波器阶数N和3dB截止频率wc、用butter求s域的频率响应的参数、用bilinear函数即利用双线性变换实现频率响应s域到z域的变换，然后用filter求滤波后信号。

6.3程序源码及注释

clear
[y,fs]=audioread('beautiful.wav');%读取音乐信号
y1=y(:,1);%取单列                            
N=length(y1);%求信号长度 
Y=fft(y1,N);%傅里叶变换
t=(0:(N-1))/fs;%时间
f=0:fs/N:fs*(N-1)/N;%频率
w=2*f/fs;%归一化
%原信号时域波形和频谱图
figure(1);
subplot(2,1,1);
plot(t,y1);
xlabel('时间t');
ylabel('幅度');
title('原信号时域波形');
subplot(2,1,2);
plot(w,abs(Y));
xlabel('w');
ylabel('幅度');
title('原信号频谱图');
%加三余弦混合噪声
A0=0.05;%噪声幅度
d1=A0*cos(2*pi*3000*t);%3kHz的余弦噪音
d2=A0*cos(2*pi*5000*t);%5kHz的余弦噪音
d3=A0*cos(2*pi*8000*t);%8kHz的余弦噪音
x1=d1+d2+d3;%噪声叠加
X1=fft(x1);%噪声傅里叶变换
y2=y1+(x1)';%加三余弦混合噪声后信号
Y1=fft(y2);%加噪信号傅里叶变换
figure(2);
subplot(3,2,1);
plot(w,abs(X1));
xlabel('w');
ylabel('幅度');
title('三余弦噪声频谱');
subplot(3,2,3);
plot(t,y2);
xlabel('时间t');
ylabel('幅度');
title('加三余弦噪声信号时域波形');
subplot(3,2,5);
plot(w,abs(Y1));
xlabel('w');
ylabel('幅度');
title('加三余弦噪声信号频谱图');
%加随机白噪声
x2=0.5*randn(N,1);%噪声
X2=fft(x2);%噪声傅里叶变换
y3=y1+x2;%加随机噪声后信号       
Y2=fft(y3);%加噪信号傅里叶变换
subplot(3,2,2);
plot(w,abs(X2));
xlabel('w');
ylabel('幅度');
title('随机白噪声频谱图');
subplot(3,2,4);
plot(t,y3);
xlabel('时间t');
ylabel('幅度');
title('加随机噪声信号时域波形');
subplot(3,2,6);
plot(w,abs(Y2));
xlabel('w');
ylabel('幅度');
title('加随机噪声信号频谱图');
%巴特沃斯滤波器去噪
wp=0.1;ws=0.16;Rp=1;Rs=15;
[n1,Wn1]=buttord(wp,ws,Rp,Rs,'s');%求低通滤波器的阶数和截止频率
[b1,a1]=butter(n1,Wn1,'s');%求s域的频率响应的参数
[b2,a2]=bilinear(b1,a1,1);%利用双线性变换实现频率响应s域到z域的变换
z1=filter(b2,a2,y2);%求滤三余弦混合噪声后的信号
z2=filter(b2,a2,y3);%求滤随机噪声后的信号
  %画滤噪后的图
m1=fft(z1);
m2=fft(z2);
figure(3);
subplot(2,2,1);
plot(t,z1);
xlabel('时间t');
ylabel('幅度');
title('滤三余弦噪声后信号时域波形');
subplot(2,2,3);
plot(w,abs(m1));
xlabel('w');
ylabel('幅度');
title('滤三余弦噪声后信号频谱');
subplot(2,2,2);
plot(t,z2);
xlabel('时间t');
ylabel('幅度');
title('滤随机噪声后信号时域波形');
subplot(2,2,4);
plot(w,abs(m2));
xlabel('w');
ylabel('幅度');
title('滤随机噪声后信号频谱');
%播放
sound(y1,fs);%原音乐
pause(5);
sound(y2,fs);%加三余弦噪声
pause(5);
sound(y3,fs);%加随机噪声
pause(5);
sound(z1,fs);%滤三余弦噪声
pause(5);
sound(z2,fs);%滤随机噪声

6.4运行结果

仿真结果分析：从图中可以看到三余弦混合噪声的频谱图以及信号加三余弦混合噪声后的时域波形和频谱，通过播放感受到了噪声对信号的影响；滤波之后对噪声的改善很大，噪声变小，原声音更加清晰，只是巴特沃斯滤波会把一部分原信号频率滤掉，声音会有点低沉。

7、语音/音乐信号的幅频滤波及相频分析

7.1题目要求

① 设计低通滤波器（可自行选择不同的截止频率），滤除原始语音/音乐信号的高频信息，观察滤波前后的幅度频谱，并比较滤波前后的音乐效果，感受高频信息对语音/音乐信号的影响；

② 设计高通滤波器（可自行选择不同的截止频率），滤除原始语音/音乐信号的低频信息，观察滤波前后的幅度频谱，并比较滤波前后的音乐效果，感受低频信息对语音/音乐信号的影响；

③ 选取两段不同的语音/音乐信号，分别将其幅度谱与相位谱交叉组合构成新的语音/音乐信号，播放比较组合后的音乐与原始音乐，感受相频信息对语音/音乐信号的影响。

7.2设计内容及方案

① 低通滤波器设计：我这里用了巴特沃斯低通滤波器，其中用buttord求低通滤波器的阶数和截止频率，用butter求s域的频率响应的参数，用bilinear即双线性变换法实现频率响应s域到z域的变换，用filter实现滤波。

② 高通滤波器设计：我这里用了巴特沃斯低通滤波器转高通，其中用buttord求低通滤波器的阶数和截止频率，用buttap创建巴特沃斯低通滤波器原型，用zp2tf将模拟低通变高通，用bilinear即双线性变换法实现频率响应s域到z域的变换，用filter实现滤波。

③ 幅度谱和相位谱交叉：用函数abs和angle分别取信号的幅度和相位，然后将其交叉创建新的信号。

7.3程序源码及注释

clear
[y,fs]=audioread('突然好想你.wav');%读取音乐信号
y1=y(:,1);%取单列                            
N=length(y1);%求信号长度
Y=fft(y1,N);%傅里叶变换
t=(0:(N-1))/fs;%时间
f=0:fs/N:fs*(N-1)/N;%频率
w=2*f/fs;%归一化
%IIR低通滤波器
wp=0.1;ws=0.2;Rp=1;Rs=15;
[n1,Wn1]=buttord(wp,ws,Rp,Rs,'s');%求低通滤波器的阶数和截止频率
[b1,a1]=butter(n1,Wn1,'s');%求s域的频率响应的参数
[b2,a2]=bilinear(b1,a1,1);%利用双线性变换实现频率响应s域到z域的变换
z1=filter(b2,a2,y1);%滤波后信号
m1=fft(z1);
figure(1);
subplot(2,1,1);
plot(w,abs(Y));
xlabel('w');
ylabel('幅度');
title('原信号频谱');
subplot(2,1,2);
plot(w,abs(m1));
xlabel('w');
ylabel('幅度');
title('低通滤波信号频谱');
%高通滤波器
wp=0.1;ws=0.2;rp=1;rs=15;
[n,Wn]=buttord(wp,ws,rp,rs,'s');%设计模拟滤波器
[z,p,k]=buttap(n);%创建巴特沃斯低通滤波器原型
[b,a]=zp2tf(z,p,k);%由零极点转换为传递函数
[b11,a11]=lp2hp(b,a,Wn);%模拟低通变高通
[b22,a22]=bilinear(b11,a11,0.5);%利用双线性变换实现频率响应S域到Z域的变换
z2=filter(b22,a22,y1);%滤波后信号
m2=fft(z2);
figure(2);
subplot(2,1,1);
plot(w,abs(Y));
xlabel('w');
ylabel('幅度');
title('原信号频谱');
subplot(2,1,2);
plot(w,abs(m2));
xlabel('w');
ylabel('幅度');
title('高通滤波信号频谱');
%幅度谱相位谱交叉
[x,fs0]=audioread('beautiful.wav');
x0=x(:,1); 
x1=x0(1:200000);%信号2截取长度200000
y2=y1(1:200000);%信号1截取长度200000
N=length(x1);
X=fft(x1,N);
Y=fft(y2,N);
t0=(0:(N-1))/fs0;
f0=0:fs0/N:fs0*(N-1)/N;
w0=2*f0/fs0;%数字角频率
%画两信号的幅度谱和相位谱
figure(3);                                 
subplot(2,2,1);
plot(w0,abs(Y));
xlabel('w');
ylabel('幅度');
title('信号1的幅度谱');
subplot(2,2,2);
plot(w0,abs(X));
xlabel('w');
ylabel('幅度');
title('信号2的幅度谱');
subplot(2,2,3);
plot(w0,angle(Y));   
xlabel('w');
ylabel('相位');
title('信号1的相位谱');
subplot(2,2,4);
plot(w0,angle(X));
xlabel('w');
ylabel('相位');
title('信号2的相位谱');
F1=abs(Y).*exp(1i*angle(X));%音乐信号1的幅度谱和音乐信号2的相位谱交叉
F2=abs(X).*exp(1i*angle(Y));%音乐信号2的幅度谱和音乐信号1的相位谱交叉y3=real(ifft(F1,N));%对交叉得到的频域信号做傅里叶逆变换
y4=real(ifft(F2,N));%对交叉得到的频域信号做傅里叶逆变换
%绘制交叉得到的信号的波形和频谱
figure(4)
subplot(2,2,1);
plot(t0,y3);
xlabel('时间');
ylabel('幅度');
title('1幅度谱和2相位谱交叉信号时域波形');
subplot(2,2,2);
plot(w0,abs(F1));
xlabel('w');
ylabel('幅度');
title('1幅度谱和2相位谱交叉信号频谱');
subplot(2,2,3);
plot(t0,y4);
xlabel('时间');
ylabel('幅度')
title('2幅度谱和1相位谱交叉信号时域波形');
subplot(2,2,4);
plot(w0,abs(F2));
xlabel('w');
ylabel('幅度');
title('2幅度谱和1相位谱交叉信号频谱');
%播放
sound(y1,fs);%原音乐
pause(5);
sound(z1,fs);%低通滤波声音
pause(5);
sound(z2,fs);%高通滤波声音
pause(5);
sound(y3,fs);%幅度谱1与相位谱2交叉
pause(5);
sound(y4,fs);%相位谱1与幅度谱2交叉

7.4运行结果

仿真结果分析：通过观察原信号的归一化频谱，确定巴特沃斯高通滤波器的参数wp,ws的值并实现滤波，从低通滤波和高通滤波的频谱图中可以看出：低通滤波器滤掉了信号的高频部分，声音变得低沉；高通滤波器滤掉了信号的低频部分，声音变得尖锐。幅度谱与相位谱交叉时，通过听交叉后的语音让我感受到了相频特性对一个信号的影响，音乐幅度谱没变相位谱变还会有原声，只是整体节奏改变。

8、语音信号变声处理系统

8.1题目要求

① 电视台经常针对某些事件的知情者进行采访，为了保护知情者，经常改变说话人的声音，请利用所学的知识，将其实现；

② 要求处理后的语音信号基本不影响正常收听与理解。

8.2设计内容及方案

变声主要通过调用voice函数来实现，改变voice函数的参数，小于1时偏女声，大于1时偏男声。

8.3程序源码及注释

 

clear
[y1,fs]=audioread('beautiful.wav');%读取音乐信号
y1=y1(:,1);%取单列
N1=length(y1);%求信号长度
n=0:N1-1;%所有元素
t=n/fs;%时间
f=0:fs/N1:fs*(N1-1)/N1;%频率
Y1=fft(y1,N1);%傅里叶变换
%画原信号图
figure(1)
subplot(2,2,1);
plot(t,y1);
xlabel('时间t');
ylabel('幅度');
title('原信号时域波形图');
subplot(2,2,2);
plot(f,abs(Y1));
xlabel('频率f');
ylabel('幅度');
title('原信号频谱图');
y2=voice(y1,0.5);  %变声
N2=length(y2);
k=0:1:N2-1;
t=k/fs;
f=k*fs/N2;
Y2=fft(y2,N2);
%画变音后信号图
subplot(2,2,3);
plot(t,y2);
xlabel('时间t');
ylabel('幅度');
title('变声信号时域波形图');
subplot(2,2,4);
plot(f,abs(Y2));
xlabel('频率f');
ylabel('幅度');
title('变声信号频谱图');
sound(y2,fs);

 

voice函数：

function Y=voice(y1,f)      %更改采样率使基频改变 f>1降低;f<1升高?
f=round(f*1000);
d=resample(y1,f,1000);      %时长整合使语音文件恢复原来时长
W=400;
Wov=W/2;
Kmax=W*2;
Wsim=Wov;
xdecim=8;
kdecim=2;
X=d';
F=f/1000;
Ss=W-Wov;
xpts=size(X,2);
ypts=round(xpts/F);
Y=zeros(1,ypts);
xfwin=(1:Wov)/(Wov+1);
ovix=(1-Wov):0;
newix=1:(W-Wov);
simix=(1:xdecim:Wsim)-Wsim;
padX=[zeros(1,Wsim),X,zeros(1,Kmax+W-Wov)];
Y(1:Wsim)=X(1:Wsim);
lastxpos=0;
km=0;
for ypos=Wsim:Ss:(ypts-W)
xpos=round(F*ypos);
kmpred=km+(xpos-lastxpos);
lastxpos=xpos;
if(kmpred<=Kmax)
   km=kmpred;
else
ysim=Y(ypos+simix);
rxy=zeros(1,Kmax+1);
rxx=zeros(1,Kmax+1);
Kmin=0;
for k=Kmin:kdecim:Kmax
xsim=padX(Wsim+xpos+k+simix);
rxx(k+1)=norm(xsim);
rxy(k+1)=(ysim*xsim');
end
Rxy=(rxx~=0).*rxy./(rxx+(rxx==0));
km=min(find(Rxy==max(Rxy))-1);
end
xabs=xpos+km;
Y(ypos+ovix)=((1-xfwin).*Y(ypos+ovix))+(xfwin.*padX(Wsim+xabs+ovix));
Y(ypos+newix)=padX(Wsim+xabs+newix);
end
end

8.4运行结果

仿真结果分析：该程序可以很好的实现对声音信号的变声处理，变声后语速基本没有改变，可以听清讲述者想传达的内容，只是音色有了变化。

小结：一次收获很多的课设，希望对你有用。