实验目的和要求

综合应用数字信号处理的理论知识进行信号的频谱分析和滤波器设计，通过理论推导得到响应结论，利用MATLAB进行计算机仿真，加深对所学知识的理论界，融会贯通所学知识。

通过本次课程设计，掌握用MATLAB对语音信号进行分析和处理的能力，并进一步掌握MATLAB 设计数字滤波器的方法。

基础设计要求：

(1)录制或采集一段语音信号，画出原始语音信号的时域波形及频谱。

(2)应用Matlab平台给语音信号叠加高斯白噪声；画出加噪信号的时域波形及频谱。

(3)确定设计性能指标，设计IIR数字滤波器进行滤波，画出滤波器的频谱。
(4)使用设计好的滤波器对加噪信号进行滤波，得到去噪信号。画出去噪信号的时域及频谱图。回放滤波前后语音，对比滤波效果。

提高设计要求：

(1)在原始语音上叠加高频噪音或低频噪音，完成滤波器设计过程。

(2)设计两种及以上滤波器对语音信号进行滤波，给出不同滤波器滤波前后的信噪比（SNR），分析得出滤波效果最好的滤波器。

(3)使用GUI/app design设计语音去噪系统界面。

(4)查阅课本或文献使用其他滤波算法完成滤波。(谱减法、卡尔曼滤波等)

实验原理

语言去噪步骤

基于MATLAB的语音信号去噪设计主要分为以下五个步骤：

（1）语音采集：可在MATLAB平台上录入一段语音信号，也可直接导入收集好的语音信号；

（2）语音分析：绘制原始语音信号时域及频谱图，分析原始语音信号频谱特征；

（3）语音加噪：对原始语音信号叠加噪声，并绘制加噪信号时域及频谱图。常见噪声包括高斯白噪声、高频噪声、低频噪声等；

（4）滤波器设计：结合原始语音信号频谱，针对不同的叠加噪声，确定滤波器设计性能指标，设计相应滤波器，并绘制滤波器的频谱图，判断是否符合设计要求；

（5）去噪信号分析：利用设计的滤波器对加噪信号进行滤波，绘制去噪信号时域及频谱图，并播放去噪语音信号，与原始语音信号进行对比分析。

谱减法的理论基础

假设带噪语音中的噪声是加性的，带噪语音谱减去估计出的噪声谱便可以得到干净语音谱。
根据加性假设，干净语音、噪声和带噪语音的关系可以写成如下所示：
$$
y(n)=x(n)+d(n)
$$

其中, $y(n)$ 为采集到的带噪语音, $x(n)$ 为干净语音, $d(n)$ 为噪声。传换至频域:
$$
\begin{aligned}
& Y(\omega)=X(\omega)+D(\omega) \
& Y(\omega)=|Y(\omega)| e^{j \Phi_y(\omega)}
\end{aligned}
$$

其中 $|Y(\omega)|$ 为幅度谱 $e^{j \Phi_y(\omega)}$ 为信号相位
则干净语音可以通过带噪语音减去噪声谱得到:
$$
\hat{X}(\omega)=[|Y(\omega)|-\hat{D}(\omega)] e^{j \Phi_y(\omega)}
$$

但在实践中, $\hat{X}(\omega)$ 很有可能是负数, 所以常会做一个半波整流:

由于对负值进行半波整流, 导致帧频谱的随机频率上出现小的、独立的峰值, 变换到时域上面, 这些峰值听起来就像帧与帧之间频率随机变化的多颤音, 也就是通常所说的 “音乐噪声” (Musical Noise)

实验方法与内容

基础设计

我们选取了一段CCTV的音频，音频的读取通过audioread函数实现
由于我们选取的音频为双通道的，为了后续处理方便将其改为单通道

y=y(:,1);

通过plot函数即可得到它的时域图，通过fft即可得到它的频谱图，通过randn函数来生成白噪声

1	x=randn(m,1); %m为CCTV音频的长度

将其与原信号相加即可得到带噪信号，再对带噪信号用plot函数即可画出时域图，通过fft得到它的频谱图

通过观察频谱图的信号和噪声分布，选择合适的性能指标，通过buttord函数得到相应的N和Wn

1	[N,Wn]=buttord(Wp,Ws,ap,as);

通过butter函数得到传递函数的系数B和A

1	[B,A]=butter(N,Wn);

用filter函数得到经降噪后的音频

1	y2=filter(B,A,y1);

然后通过sound函数播放降噪前后的音频，比对得到滤波效果的好坏并调整性能指标

提高设计

非GUI设计

高频或低频噪声的获取可通过对高斯白噪声通过一个高通或低通滤波器得到，由于我们选择的音频信息都集中在低频段，若添加低频噪声则用传统的滤波器很难消除，所以我们选择在原信号上添加高频噪声

1 2	[N,wc]=buttord(0.8,0.5,3,35); [B,A]=butter(N,wc,'high');

高斯白噪声通过这样一个高通滤波器即可变为高频信号，由于信号在低频，噪声在高频，故低通滤波器能够很好的满足我们的需求，而通过基础设计我们得到了一个巴特沃斯低通滤波器，则可以直接沿用该滤波器进行滤波

我们还可以分别设计切比雪夫Ⅰ型滤波器，切比雪夫Ⅱ型滤波器，和椭圆滤波器来完成滤波

切比雪夫Ⅰ型滤波器可通过cheb1ord函数和cheby1函数得到

[N,Wn]=cheb1ord(Wp,Ws,ap,as);
[B,A]=cheby1(N,ap,Wn);
切比雪夫Ⅱ型滤波器可通过cheb2ord函数和cheby2函数得到
[N,Wn]=cheb2ord(wp,ws,ap,as);
[B,A]=cheby2(N,Rs,Wn);

椭圆滤波器可通过ellipord函数和ellip函数得到

1 2	[N,Wn]=ellipord(wp,ws,Rp,Rs); [B,A]=ellip(N,Rp,Rs,Wn);

为了比较滤波器的滤波效果，我们需要计算滤波前后的信噪比，而我们可以选取的带噪信号相同，有相同的信噪比，只需计算滤波后的信噪比大小就能分析滤波器的性能好坏

使用滤波器进行滤波后的信噪比计算方法为：

1	SNR=20*log10(norm(y_flitered-x_flitered)/norm(x_flitered)) %y_flitered为滤波后的带噪信号,x_flitered为滤波后的噪声,norm为2范数

即滤波后的带噪信号减去滤波后的噪声即为经滤波后的纯净信号

GUI设计

我们在界面中放入三个坐标轴，分别用于观察原信号，加噪信号去噪信号的时域和频域图。

放入两个列表框，用于选择添加信号的类型和滤波器的类型或者滤波方法。

第一个列表框的选项如下

第二个列表框选项如下

列表框的选择选项功能用switch语句实现

一个可编辑文本框用于设置加噪信号的信噪比，输入即为SNR

1 2	SNR=0; %信噪比大小 x=x/norm(x,2).10^(-SNR/20)norm(y); %x为噪声,y为原信号

四个可编辑文本框用于输入滤波器的性能指标

选用传统滤波器滤波后会计算出相应的信噪比

若选用谱减法进行滤波，由于该方法滤波后的带噪信号减去滤波后的噪声不为纯净信号，则无法通过该方法计算信噪比

经计算得到的信噪比输出在可编辑文本框中

有三个播放按钮用于播放原信号，加噪信号，降噪信号，通过比对分析滤波效果的好坏

实验原始纪录

[y,Fs]=audioread("C:\Users\kiana\Desktop\cctv.wav");
% y=y(:,1);
m=length(y);
x=randn(m,1);
SNR=0;                 %信噪比大小
x=x/norm(x,2).*10^(-SNR/20)*norm(y);
y1=x+y;
% sound(y,Fs);
subplot(4,2,1);plot(y);
xlabel('时间');
ylabel('幅值');
title('原始语音信号');
subplot(4,2,3);plot(y1)
xlabel('时间');
ylabel('幅值');
title('加噪语音信号');
% sound(y1,Fs);
% audiowrite('C:\Users\kiana\Desktop\cctv_noise.wav',y1,Fs); 
Y = fft(y);  % 傅里叶变换
Y1=fft(y1);
P = abs(Y);  % 双边频谱
P1=abs(Y1);
n=0:m-1;
w=2*n/m;
subplot(4,2,2);plot(w,P);title('原始语音频谱');
subplot(4,2,4);plot(w,P1);
xlabel('频率\omega/\pi');
ylabel('幅度');
title('加噪语音频谱');
fp=0.15;fs=0.24;ap=3;as=35;
[N,Wn]=buttord(fp,fs,ap,as);
[B,A]=butter(N,Wn);
% [N,Wn]=cheb1ord(fp,fs,ap,as); %使用其他滤波器即注释当前的滤波器和去相应滤波器的注释
% [B,A]=cheby1(N,ap,Wn);
% [N,Wn]=cheb2ord(fp,fs,ap,as);
% [B,A]=cheby2(N,ap,Wn);
[H,W]=freqz(B,A);
% plot(abs(H));
% [N,Wn]=ellipord(fp,fs,ap,as);
% [B,A]=ellip(N,ap,as,Wn);
y2=filter(B,A,y1);
x2=filter(B,A,x);
Ps=sum((y2-x2).^2);
Pu=sum(x2.^2);
SNR1=10*log10(Ps/Pu); 
subplot( 4,2,5);plot(y2);xlabel('时间');
ylabel('幅值');title("去噪语音信号");
Y2=fft(y2);
P2=abs(Y2);
subplot(4,2,6);plot(w,P2);xlabel('频率\omega/\pi');
ylabel('幅度');title('去噪语音频谱');
subplot(4,2,7);plot(W/pi,20*log10(abs(H)));
% axis([0,1,-440,0]);title("巴特沃斯低通滤波器");
xlabel("\omega/pi");ylabel("增益(dB)");
% sound(y2,Fs);
% audiowrite('C:\Users\kiana\Desktop\cctv_denosied.wav',y2,Fs);

实验结果及分析

基础部分

我们选取的音频的时域和频域图如下

可以看到它的信号主要集中在低频部分

添加高斯白噪声得到信噪比为0的带噪信号，它的时域和频域图如下

发现原来各频率上均已有值，噪声的幅度均匀分布在各频率段，通过观察我们发现，在$0.15$之前的频率就包含了原信号频谱的绝大部分，而在$0.24$之后几乎没有值，所以我们可以设计相应的滤波器以保留低频段信息，滤去其他信息

得到的滤波器响应曲线

经过滤后的信号的时域和频域图如下

通过先前的信噪比计算方法得到滤波后的信号信噪比约为$7.87$这是一个合理的值，由于加噪信号信噪比为 $0$，即噪声功率$Pu$和纯净信号功率$Ps$相等，我们通过计算可得它的值应该大于$10\times log10(Ps/0.24Pu)=6.19$,小于$10\times log10（Ps/0.15Pu）=8.24$