基于MATLAB的语音信号分析与处理

摘要

        数字滤波技术是数字信号处理的一个重要组成部分, 滤波器的设计是信号处理的核心问题之一。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。本篇设计报告所选主题是基于MATLAB的IIR数字滤波器的设计与实现。MATLAB功能强大、简单易学、编程效率高,深受广大科技工作者的欢迎。MATLAB具有信号分析工具箱,不需具备很强的编程能力,就可以很方便地进行信号分析、处理和设计，利用MATLAB信号处理工具箱可以快速有效地设计各种数字滤波器。本次课程设计需要综合运用数字信号处理的理论知识对加噪声语音信号进行时域、频域分析和滤波，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现。在设计实现的过程中，使用巴特沃斯、切比雪夫原型低通滤波器通过双线性变法设计IIR数字低通、带阻滤波器，并利用MATLAB函数作为辅助工具完成设计中的计算与图形的绘制。通过对对所设计滤波器的仿真和频率特性分析，可知利用MATLAB信号处理工具箱可以有效快捷地设计IIR数字滤波器，过程简单方便。最后在语音信号噪声滤除方面的效果非常明显，噪声得到有效滤除。

关键词：语音信号处理;IIR滤波器设计;GUI设计;双线性变换法;巴特沃斯滤波器

3.1 语音信号的采集及频谱分析

利用录音软件用单声道录制一段自己的声音，时间在2S内。然后MATLAB软件平台下，利用函数audioread对语音信号进行采样，记住采样频率和采样点数。为了将原始的模拟语音信号转变为数字信号，必须进行采样和量化，进而得到时间和幅度上均为离散的数字语音信号。又因为 MATLAB只能处理“.wav”结尾的音频信号，所以在录制声音的时候，格式不匹配就需要利用格式工厂软件进行格式转化。

将采集的信号进行频谱分析，以下是其时域及其频域波形，如图3-1所示。从图中可以看出语音信号频率一般在300-3000Hz左右。而我在程序中有一步是y=y(:,1);目的是取语音信号y的第一列,考虑到可能处理的是立体声音频文件，所以先进行这个操作是必要的。

3.2 语音信号加噪处理

本次实验添加了两种高频正弦噪声和两种高斯白噪声，共添加了四个噪声，分别是10800Hz,10800Hz+14800Hz,高斯白噪声，带限高斯白噪声。

在MATLAB中产生高斯白噪声非常方便，我们可以直接应用两个函数：一个是WGN，另一个是AWGN。WGN用于产生高斯白噪声，AWGN则用于在某一信号中加入高斯白噪声。也可直接用randn函数产生高斯分布序列。在本次设计中，我们是利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。Randn函数有两种基本调用格式：Randn(n)和Randn（m,n）,前者产生n×n服从标准高斯分布的随机数矩阵，后者产生m×n的随机数矩阵。在这里，我们选用Randn（m,n）函数。但是考虑到高斯白噪声无法全部滤除，所以又将高斯白噪声经过一个高通滤波器得到一个其高频部分，以免在滤除时效果不佳。

而正弦噪声则可直接在初始语音信号上加上一个正弦函数如noise=0.1*sin(2*pi*f*k/fs)；当然我们也可以添加多个不同频率的正弦噪声这里以两个为例，noise=0.1*sin(2*pi*14800*k/fs)+0.1*sin(2*pi*10800*k/fs);

以上四种加噪信号的频谱图如图3-2，3-3，3-4，3-5所示。在图3-2中可以发现，在频率为10800Hz处即为添加的噪声频谱。在图3-3中，在频率为10800Hz处和14800Hz处时添加的噪声的频谱。在图3-4中可以发现，整段频谱都加上了噪声，因此无法全部滤除，效果不佳，而在图3-5中，在9000Hz开始加上高斯噪声，通过低通滤波器可以将其尽可能的滤除。因为我们初始的语音信号为列向量，因此需要对噪声进行转置才能加在语音信号上，即noise=noise';

带限高斯白噪声实则就是将一个高斯白噪声经过一个高通滤波器后得到的，目的是不影响原语音信号的信息频谱，下图为带限高斯白噪声的时域波形和频谱，如图3-6。

3.3设计IIR数字滤波器

我使用双线性变换法设计上面的三个滤波器。在 Matlab中，可以利用函数butter， cheby1和ellip设计IIR滤波器；利用Matlab中的函数freqz画出各滤波器的频率响应。我使用的是巴特沃斯和切比雪夫模拟滤波器，通过双线性变换法设计数字滤波器，三种滤波器的幅频特性如图3-6，图3-7，图3-8。

无论是数字滤波器还是模拟滤波器，他们技术指标的建立都是以所谓的“固有衰减”参数为参照。以数字滤波器为例，固有衰减参数（关于之旅分量归一化）定义为

                    A(w)=-20lg|H(ejw)|-20lg|H（ej0）|            （式3-1）

目前IIR数字滤波器设计的最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法，它不但有完整的设计公式，而且还有较为完整的图表供查询，因此，充分利用这些已有的资源将会给数字滤波器的设计带来很大方便。

3.3.1双线性变换法

脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。

对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：

（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。

（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为w=2/T tan(0.5Ω)。

（3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。

（4）设计模拟低通滤波器。

（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。

（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。

3.3.2性能指标

下面是实验所用滤波器的性能指标：

1）低通滤波器1性能指标:fp=8500Hz,fs=10500Hz,Ap=1dB,As=50dB。

2) 低通滤波器2性能指标:fp=2000HZ,fc=3500HZ,As=50dB,Ap=1dB。

3）带阻滤波器性能指标：fp1=10500Hz,fp2=11100Hz，fs1=10750Hz，fs2=10850Hz，Ap=1dB,As=50dB。

我使用的是巴特沃斯和切比雪夫模拟滤波器，通过双线性变换法设计数字滤波器，三种滤波器的幅频特性如图3-6，图3-7，图3-8。

3.4 滤波去噪及去噪信号频谱分析

 

用上述设计的滤波器分别对加不同噪声的语音信号进行滤波，并对语音信号加入噪声后的频谱分析。在MATLAB中，IIR滤波器利用filter对信号进行滤波。 比较滤波前后语音信号的波形及频谱，对滤波器性能进行频谱分析。在一个窗口先后画出滤波前后语音信号的波形及频谱。

图3-9为加单频正弦噪声经过带阻滤波器的去噪信号时域波形，我们可以发现，频率在10700-10900Hz范围内的信号被滤除。经过加噪处理后的语音信号，可在Matlab中用函数sound对声音进行回放。其调用格式：sound(y,Fs)，sound(y)和sound(y,Fs,bits)。可以察觉滤波前后的声音有明显的变化。通过对比前后初始语音信号和去噪语音信号，发现对噪声去噪很完整，去噪效果很好，当然具体指标得通过计算信噪比得到。

图3-10为加两个正弦噪声（10800Hz和14800Hz）的语音信号经过带阻滤波器的去噪信号时域波形，我们可以发现，频率在高于8500Hz之外的信号被滤除，所加的10800Hz和14800Hz噪声也被很好的滤除。通过对比前后初始语音信号和去噪语音信号，发现对噪声去噪很完整，去噪效果很好，具体指标得通过计算信噪比得到。

图3-11为加高斯白噪声的语音信号经过低通滤波器的去噪信号时域波形，我们可以发现，频率在高于2500Hz之外的信号几乎被滤除，所加的2500Hz以外的高斯白噪声也被很好的滤除。但是由于给整段频谱都加上了高斯白噪声，因此不可能全部滤除，仅能去除大部分噪声，改善信号质量。具体去噪效果指标可以通过计算信噪比得到。同时，由于部分语音信号高频分量被滤除，去噪信号听起来较闷，声音较为低沉。

由于加整段高斯白噪声无法全部滤除，我将一个高斯白噪声通过一个高通滤波器，留下频率高于9500Hz的噪声，再到语音信号上，使其不会影响低频的语音信号。图3-12为加带限高斯白噪声的语音信号经过低通滤波器的去噪信号时域波形。我们可以发现，频率在高于9500Hz之外的信号几乎被滤除，所加的高斯白噪声也被很好的滤除。由于给频谱高频部分加上了高斯白噪声，因此可以通过低通滤波器全部滤除，有效改善信号质量。具体去噪效果指标可以通过计算信噪比得到。

3.5设计GUI界面

为了使编写的程序操作方便，设计用户界面。在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。

3.5.1系统界面设计工具—GUI概述

图形用户界面（graphical user interfaces ，GUI）则是由窗口、光标、按键、菜单、文字说明等对象（objects）构成的一个用户界面。用户通过一定的方法（如鼠标或键盘）选择、激活这些图形对象，使计算机产生某种动作或变化，比如实现计算、绘图等。在MATLAB中GUI是一中包含多种对象的图形窗口，并为GUI开发提供一个方便高效的集成开发环境GUIDE。GUIDE主要是一个界面设计工具集，MAYLAB将所有GUI支持度控件都集成在这个环境中，并提供界面外观、属性和行为响应方式的设置方法。GUIDE将设计好的GUI保存在一个FIG文件中，同时生成M文件框架。

FIG文件：包括GUI图形窗口及其所有后裔的完全描述，包括所有对象属性的属性值。它是一个二进制文件调用hsave课保存图形窗口时将生车该文件。M文件包括GUI设计、控件函数以及定义为子函数的用户控件回调函数，主要用于控制GUI展开时的各种特征。 GUI创建包括界面设计和控件编程两部分，主要步骤如下。第一步：通过设置GUIDE应用程序的选项来运行GUIDE；第二步：使用界面设计编辑器进行面设计；第三步：编写控件行为响应控制（即回调函数）代码。

3.5.2创建GUI界面

首先我们新建一个GUI文件：File/New/GUI , 选择Blank GUI(Default)。

3.5.3创建开关按钮

创建8个静态文本（Static text），18个Push Button，3个UIpanel如下图所示（属性可双击更改）

参考文献

[1]高西全、丁玉美.数字信号处理(第四版)[M].西安：西安电子科技大学出版社，2016

[2]汪依帆,王敏,李海龙.基于 MATLAB 的 IIR 数字滤波器设计研究[J].数字通信世界,2019(8):104.

[3]文胜福、王硕、魏宇豪、王浩.基于IIR数字滤波器的人体脉搏信号分析[A].重庆：重庆邮电大学，2020

[4]潘一飞，王家琪，陈波. 电力系统谐波检测中数字滤波器的设计[A].河北：华北理工大学，2020

[5]ANSI Sl.42 -2001 Design response of weighting networks foracoustical measurements[S]. New York：Acoustical Society of America，2001.

[6]唐求，吴娟，邱伟，沈洁，滕召胜. 基于 MICA 的声级计频率计权数字 IIR 滤波器设计[A].湖南：湖南大学，2020

[7]Sheng-Ming Yang, Shuenn-Jenn Ke.Performance Evaluation of a Veloci ty Observer for Accurate Velocity Esti mation of Servo Motor Drives[J].IEEE Trans.On Industrial Aplications, 2000, 36 (1) :98-104

[8]杜仁慧,陶春荣,张伟,施永柱.IIR数字滤波器在电机转速测量中的应用[J].雷达与对抗,2017,37(02):64-68.

[9]邢子龙,闫利文,杨振峰,赵小超.IIR数字滤波器在砌块成型机振动信号滤波中的应用[J].现代机械,2015(06):53-55.

[10]林竞连. IIR滤波器的分析及音频提取[A].西北：西北民族大学，2021

%hObject：这个参数是一个图形界面元素的句柄，它指向被用户操作的那个图形界面元素。
% 例如，如果用户点击了一个按钮，那么回调函数中的hObject参数就会指向这个按钮。
%eventdata：这个参数包含关于事件的数据信息。具体的内容取决于发生的事件。
% 例如，如果用户点击了一个按钮，eventdata可能包含有关点击的信息，如点击的位置和时间。
%handles：这个参数是一个包含所有图形界面元素句柄的结构体。通过这个参数，我们可以访问所有的图形界面元素。
% 在回调函数中，我们可以通过handles参数来获取其他图形界面元素的信息，并根据需要更新它们的状态。
function varargout = wxk(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @kesegui_OpeningFcn, ...
                   'gui_OutputFcn',  @kesegui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
function kesegui_OpeningFcn(hObject, ~, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
function varargout = kesegui_OutputFcn(~, ~, handles) 
varargout{1} = handles.output;
function edit1_CreateFcn(hObject, ~, ~)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit2_CreateFcn(hObject, ~, ~)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function uibuttongroup4_SelectionChangedFcn(~, ~, ~)
% --- Executes on button press in radiobutton14.
function radiobutton14_Callback(~, ~, ~)
%********************************************************************** %
function axes1_CreateFcn(~, ~, ~)
function axes1_ButtonDownFcn(~, ~, ~)
function uibuttongroup3_CreateFcn(~, ~, ~)
function axes2_CreateFcn(~, ~, ~)
function axes2_ButtonDownFcn(~, ~, ~)
function pushbutton8_ButtonDownFcn(~, ~, ~)
function pushbutton16_ButtonDownFcn(~, ~, ~)
function pushbutton16_KeyPressFcn(~, ~, ~)
function pushbutton14_KeyPressFcn(~, ~, ~)
function pushbutton14_ButtonDownFcn(~, ~, ~)
function pushbutton14_CreateFcn(~, ~, ~)
function pushbutton15_KeyPressFcn(~, ~, ~)
function pushbutton15_ButtonDownFcn(~, ~, ~)
function pushbutton15_CreateFcn(~, ~, ~)
function pushbutton19_KeyPressFcn(~, ~, ~)
function pushbutton19_ButtonDownFcn(~, ~, ~)
function pushbutton19_CreateFcn(~, ~, ~)
function pushbutton25_DeleteFcn(~, ~, ~)
%****************************************************************%
% 录音
function pushbutton17_Callback(~, ~, ~)
Fs=48000;                  %采样频率
duration=4;                %时间长度
nBits = 16 ;               %采样位数
nChannels = 2 ;            %通道数
ID = -1; %默认音频输入设备
recObj = audiorecorder(Fs,nBits,nChannels,ID);
disp('Start Recording!')
recordblocking(recObj,duration);
disp('End of Recording!');
%获取录音数据
Recording = getaudiodata(recObj);
%保存录音数据
filename = './2.wav'; 
audiowrite(filename,Recording,Fs);

%导入音频
function pushbutton1_Callback(~, ~, handles)
global y fs  %定义全局变量
[filename,pathname]=uigetfile({'*.wav;*.m4a;*.mp3;*.mov'},'选择语音信号');
if isequal(filename,0)||isequal(pathname,0)
    return;
end
yy=[pathname,filename];%从文件读取音频数据
[y,fs]=audioread(yy);
axes(handles.axes1);

%播放
function pushbutton5_Callback(~, ~, ~)
global y fs %定义全局变量
sound(y,fs);

%原语音时域图
function pushbutton7_Callback(~, ~, handles)
global y N  t   %定义全局变量
N=length(y);%求出采样值的长度
t=0:N-1;%时域范围
axes(handles.axes1);
plot(t,y);
title('原语音信号的时域图');
xlabel('时间');ylabel('振幅');

%原语音时域图频域图
function pushbutton16_Callback(~, ~, handles)
global y fs N  k D  %定义全局变量
Y=fft(y,N); %做傅里叶变换
k=0:N-1;
D=fs/N;%计算谱线间隔
axes(handles.axes1);
stem(k*D,abs(Y),'Marker','none');axis([0,5000,0,8000]);
title('原语音信号的频谱');
xlabel('频率');ylabel('幅值');

% 添加正弦噪声
function pushbutton14_Callback(~, ~, ~)
global  fs  noise  k %定义全局变量
f=10800; %噪声频率
noise=0.1*sin(2*pi*f*k/fs); %设定要加进去的噪声0.1*sin(2*pi*f)
noise=noise';%转置维度行列y为列向量，所以noise要转置

%添加两个正弦噪声
function pushbutton15_Callback(~, ~, ~)
global  fs  noise  k %定义全局变量
noise=0.1*sin(2*pi*14800*k/fs)+0.1*sin(2*pi*10800*k/fs); %设定要加进去的噪声0.01*sin(2*pi*f)
noise=noise';%转置  维度  行 列 

%加噪后时域图
function pushbutton9_Callback(~, ~, handles)
global y fs N noise xa  t  %定义全局变量
t=0:N-1;
xa=y+noise; %加噪声
axes(handles.axes2);
plot(t,xa);
title('加噪声后信号的时域图');
sound(xa,fs);%放声音

%加噪后频域图
function pushbutton8_Callback(~, ~, handles)
global y fs N noise xa Xa  k D  %定义全局变量
xa=y+noise; %加噪声
k=0:N-1;
D=fs/N;%计算频率分辨率,计算谱线间隔
Xa=fft(xa,N); %做傅里叶变换
axes(handles.axes2);
stem(k*D,abs(Xa),'Marker','none'); axis([0,15000,0,8000]);
title('加噪声后的语音信号的频谱')%在f=1800hz处有了噪声信号


%巴特沃斯带阻去除正弦噪声（IIR）
function pushbutton18_Callback(~, ~, handles)
global  fs N  xa   yx Y %定义全局变量
T=2;
fp1=10500;fp2=11100;
fs1=10750;fs2=10850;
wp1=2*pi*fp1/fs;wp11 = 2/T*tan(wp1/2);%前者为数字角频率，后者为模拟带阻角频率
wp2=2*pi*fp2/fs;wp22 = 2/T*tan(wp2/2);
ws1=2*pi*fs1/fs;ws11 = 2/T*tan(ws1/2);
ws2=2*pi*fs2/fs;ws22 = 2/T*tan(ws2/2);
wp=[wp11, wp22];
ws=[ws11 ,ws22];
[n13,wn13]=cheb1ord(wp,ws,1,50,'s');%模拟带阻角频率得到原型低通（切比雪夫）阶数和截频 
[b13,a13]=cheby1(n13,1,wn13,'stop', 's');%原型低通转换成模拟带阻滤波器（系统函数）
[num11, den11] = bilinear(b13, a13,0.5); %模拟带阻由双线性变换得到数字带阻滤波器 FS=1/T，T=2
[H,w]=freqz(num11, den11);
dbH=20*log10(abs(H));
axes(handles.axes3);%句柄，子图显示
plot(w/2/pi*fs,dbH,'LineWidth',2);
title('带阻的幅频特性图');
axis([0,12000,-40,2]);
yx=filter(num11,den11,xa);
Y=fft(yx,N);%频域相乘

% 低通-平滑
function pushbutton19_Callback(~, ~, handles)
global  fs N  xa  b a yx  Y %定义全局变量
Fp =8500; 
Fs =10500; 
T=2;
wp = 2*pi*Fp/fs;
ws = 2*pi*Fs/fs;
wp1 = 2/T*tan(wp/2);
ws1 = 2/T*tan(ws/2);
[n, wc] = buttord(wp1, ws1, 1, 50, 's'); 
[b, a] = butter(n, wc, 's'); 
[num11, den11] = bilinear(b, a,0.5); % FS=1/T，T=2
[h,w]=freqz(num11, den11);
h1 = 20*log10(abs(h));
axes(handles.axes3);%句柄
plot(w/2/pi*fs,h1,'LineWidth',2);title('低通的幅频特性图') ;
yx=filter(num11, den11,xa);
Y=fft(yx,N);%频域相乘

%去噪后时域
function pushbutton20_Callback(~, ~, handles)
global  fs  t  yx  %定义全局变量,yx为滤波后信号
axes(handles.axes3);%句柄
sound(yx,fs);
%t=0:size(yx)-1;
plot(t,yx);

%去噪后频域
function pushbutton21_Callback(~, ~, handles)
global  N  k D  yx Y %定义全局变量
Y=fft(yx,N);%傅里叶变换查看频率
axes(handles.axes3);%子图3
stem(k*D,abs(Y),'Marker','none');title('去噪后');
axis([0,15000,0,8000]);

%退出
function pushbutton4_Callback(~, ~, ~)
clc
close all
close(gcf)

%高斯白噪声
function pushbutton26_Callback(~, ~, ~)
global y  noise  %定义全局变量
y=y(:,1); %第一列
y=y';
noise=0.7*rand(1,length(y));

% 小低通
function pushbutton27_Callback(~, ~, handles)
global  fs N  xa  b a yx  Y %定义全局变量
 Fp =2000; 
 Fs =3500; 
 T=2;
 wp = 2*pi*Fp/fs;
 ws = 2*pi*Fs/fs;
 wp1 = 2/T*tan(wp/2);
 ws1 = 2/T*tan(ws/2);
 [n, wc] = buttord(wp1, ws1, 1, 50, 's'); 
 [b, a] = butter(n, wc, 's'); 
 [num11, den11] = bilinear(b, a,0.5); % FS=1/T，T=2
 [h,w]=freqz(num11, den11);
 h1 = 20*log10(abs(h));
 axes(handles.axes3);%句柄
 plot(w/2/pi*fs,h1,'LineWidth',2);title('低通的幅频特性图') ;
 yx=filter(num11, den11,xa);
 Y=fft(yx,N);%频域相乘

%停止播放
% --- Executes on button press in pushbutton28.
function pushbutton28_Callback(~, ~, ~)
clear sound

%窄带高斯白
function pushbutton29_Callback(~, ~, ~)
global y  noise  %定义全局变量
y=y(:,1); %第一列
y=y';
noise1 = 0.5*randn(size(y)); % 生成与原始信号大小一致的高斯白噪声信号   
Ft=20000;
Fp=100000;    
Fs=120000;    
wp1=tan(pi*Fp/Ft);    
ws1=tan(pi*Fs/Ft);        
wp=1;    
ws=wp1*wp/ws1;    
[n13,wn13]=cheb1ord(wp,ws,1,50,'s');    
[b13,a13]=cheby1(n13,1,wn13,'s');       
% 确保截止频率是非负的
wn13_abs = abs(wn13);
[num,den]=lp2hp(b13,a13,wn13_abs);
[num13,den13]=bilinear(num,den,0.5);
if num(1) > 0 && den(1) > 0  
    noise= filter(num13,den13,noise1);                                                     
else  
    error('Filter is not stable');  
end   

function pushbutton30_Callback(~, ~, ~)
global  fs  noise  k %定义全局变量
f=10800; %噪声频率
noise=0.1*sin(2*pi*f*k/fs); %设定要加进去的噪声0.1*sin(2*pi*f)
noise=noise';%转置维度行列