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('回声时域波形');
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);%播放低频率音乐
3 语音/音乐信号的抽取（减抽样）
3.1 题目要求
① 观察语音/音乐信号频率的上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；
② 输出减抽样语音/音乐信号的波形和频谱，观察现象，给出理论解释；
③ 播放减抽样语音/音乐信号，注意抽样率改变，比较不同抽取间隔下的声音，解释现象。
3.2 设计内容及方案
本题研究的唯一内容就是对信号以不同间隔抽样，这里抽取了200000长度的信号，便于以一定间隔抽样。
首先以小间隔2抽样，信号声音基本和原声差不多，没有发生混叠；而后又用大间隔20抽样，信号声音有了很明显的变化，即发生了混叠。
原因：采样频率fs必须大于等于2fc，采样频率小于2fc时，发生混叠，若采样频率越小，则混叠越明显。音频fc不到5000Hz,间隔2采样频率fs为22000Hz大于两倍的fc所以不混叠；间隔20采样频率fs为2200Hz小于两倍的fc所以发生混叠。
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为间隔减抽样信号播放(混叠)
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);%高频调制播放
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);%布莱克曼窗滤波后的声音
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);%滤随机噪声
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交叉
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
