语音数字信号处理系统设计（含matlab程序）

1 概述
- 1.1.设计目的
- 1.2.设计要求
- 1.3.功能描述
- - 1.3.1数字语音信号采集
  - 1.3.2时域分析与频谱分析
  - 1.3.3叠加噪声
  - 1.3.4小波去噪
  - 1.3.5数字滤波器处理
  - 1.3.6设计框图
2 原理分析
- 2.1.语音的产生和感知
- 2.2.语音信号的数学模型
- - 2.2.1激励模型
  - 2.2.2声道模型
  - 2.2.3辐射模型
- 2.3.语音信号的时域、频域特性和短时分析技术
- - 2.3.1语音信号分帧
  - 2.3.2短时时域分析
  - 2.3.3短时频域分析
- 2.4.带噪语音处理
- - 2.4.1带噪语音的产生
  - 2.4.2小波变换去噪
  - 2.4.3 IIR数字滤波器设计
3 matlab程序
4 参考文献

1 概述

主要通过MATLAB GUI建立了一个简易的语音数字信号处理系统，可以实现数字语音信号采集、时域分析、频谱分析、加噪、降噪等功能，同时利用双线性变换法设计了一个IIR滤波器，将非带限的模拟滤波器映射为最高频率为的带限模拟滤波器，然后再将模拟滤波器转换为数字滤波器，对比滤波前后的时域及频谱，并绘制出频率响应曲线。由于噪声无处不在，所以采集的原始语音信号不可避免地带有噪声，为了衡量降噪效果，需要对其先加一个已知信噪比的噪声，这里采用的是添加高斯白噪声。为了便于对语音信号的时域-频域分析需要先对原始语音信号进行快速傅里叶变换（fast Fourier transform），即FFT变换然后再绘制出时域图与频谱图作进一步分析。

1.1.设计目的

数字滤波器设计是“数字信号处理”课程的主要内容之一；由于其涉及到模拟滤波器、数字信号频域变换、周期信号的傅里叶变换、频率响应特性、抽样定理及频率混叠等理论知识，这些知识难点光学习理论是远远不够的，它必须与实践相结合，通过matlab编程实践才能更好地掌握理论知识，同时更深刻的体会数字信号处理相关理论知识的应用性特色。
所以，设计这个语音信号处理系统和数字滤波器的目的主要有：
（1）学习并掌握语音数字信号的采集与处理；
（2）学习数字滤波器的设计方法；
（3）学习语音信号的短时时域处理与短时频域处理技术；
（4）学习分析时域图与频谱图；
（5）学习数字滤波器与模拟滤波器的转换。

1.2.设计要求

设计一个简单的音频或语音信号处理系统：
（1）音乐信号中的噪声消除，方法不限；
（2）根据语音信号的特点，设计相应的数字滤波器。

1.3.功能描述

通过MATLAB GUI建立了一个简易的语音数字信号处理系统可视化界面，主要由一个控制面板与一个显示面板构成，控制面板主要是通过按钮来实现对语音数字信号的处理[1]：数字语音信号采集、时域分析、频谱分析、加噪、降噪与滤波处理；显示面板主要用来显示对语音数字信号处理前后的各种对比图。如图1-1所示。
数字语音信号采集[2][3]：通过matlab的audiorecorder函数来访问计算机上的麦克等语音设备来记录语音，并保存为wav格式的语音文件。
时域分析与频谱分析[4][5]：并对信号做1024点FFT变换[6]，用matlab绘制出语音信号的时域波形图、频率响应图与频谱图。
叠加噪声：在matlab软件平台下给原始语音信号叠加上高斯随机噪声[7][8]，并绘制加噪前后的时域图与频谱图。
消除噪声：主要通过小波函数来对加入噪声的语音信号进行小波分解[9][10]，估计噪声方差，从而得到去噪后的语音信号，并对比去噪前后的时域图。
滤波：用双线性变换法[11]设计了巴特沃斯数字低通IIR滤波器[12]对加入高斯噪声的语音信号进行滤波处理，并绘制出巴特沃斯数字低通滤波器[13][14]的频率响应曲线，及滤波前后的时域图与频谱图。
界面如图1-1所示：

上图在采集语音信号时，是利用计算机自带的麦克风，点击录音后用男声读“蓝天”“白云”，记录5秒后结束数据采集，然后绘制出原始语音信号的时域图。

1.3.2时域分析与频谱分析

信号的时域表示形式是现实世界信号的真实形式，而信号的频域分析采用傅里叶变换实现，傅里叶变换的快速算法即快速傅里叶变换[15]，所以信号的频域形式实际上是一种数学形式，但却有具体的含义和重要的用途。有些信号的时域波形很难显示其特征，比如从男生和女生发相同语音时的时域波形很难分析出其特征，以及很难从含有噪声的语音信号的时域波形中识别和去除噪声，但是在频域分析中可以很容易进行特征分析。如图1-3所示，应用 MATLAB 软件绘制出刚才记录的一段wav格式的录音的时域波形和频谱图，由频谱图可知，这段录音的频谱处于0到1000Hz频段，并且主要集中在低频段。

1.3.4小波去噪

从信号学的角度看，小波去噪是一个信号滤波的问题，而且尽管在很大程度上小波去噪可以看成是低通滤波[18]，但是由于在去噪后还能成功地保留信号特征，所以在这一点上又优于传统的低通滤波器。由此可见，小波去噪实际上是特征提取和低通滤波功能的综合，如图1-5所示，原始信号与小波去噪后的时域图比较接近，结合图1-4加入高斯白噪声后的时域图作对比，小波去噪的效果还是很明显的。

IIR数字滤波器都是由模拟滤波器的原型变换过来的，模拟滤波器的原型不仅仅有巴特沃斯滤波器，还有切比雪夫Ⅰ型滤波器[23]、切比雪夫Ⅱ型滤波器[24][25]和椭圆滤波器[26][27]。如图1-7是巴特沃斯滤波器、切比雪夫Ⅱ型滤波器和椭圆滤波器在相同指标下的频率响应曲线：

2 原理分析

从语音的产生与感知开始，到语音信号数学模型的建立，为语音信号的处理奠定了基础。根据语音信号的特点，介绍了几种常用的短时时域分析与短时频域分析技术，如短时能量与短时平均幅度、短时平均过零率和傅里叶变换等。同时还介绍了对带噪语音的处理，如小波变换，并详细介绍了IIR数字滤波器的设计。

2.1.语音的产生和感知

语音是由发声器官在大脑的控制下的生理运动产生，发音器官包括肺、气管、喉(包括声带)、咽、鼻和口等，这些器官共同形成一条形状复杂的管道，其中喉以上的部分为声道，它随发出声音的不同形状而变化；喉的部分称为声门。发声器官中，肺和器官是整个系统的能源，喉是主要的声音产生机构，而声道则对生成的声音进行调制[30]。
产生语音的能量，来源于正常呼吸时肺部呼出的稳定气流，喉部的声带既是阀门又是震动部件。二两声带间的部位为声门。说话时，声门处气流冲击声带产生震动，然后通过声道响应变成语音。发不同音时声道形状不同，所以能听到不同的声音。喉部的声带对发音影响很大，其为语音提供主要的激励源：声带震动产生声音。声带开启和闭合使得气流形成一系列脉冲。没开启和闭合一次的时间即震动周期，称为基音周期，其倒数为基因频率，简称基频[31]。
语音由声带振动或不经声带振动而产生，其中由声带振动产生的称为浊音，二不由声带振动产生的称为清音。浊音包括所有原因和一些辅音，清音包括另一部分辅音。对于浊音、清音和爆破音，其激励源不同。浊音是位于声门处的准周期麦种序列，清音是位于声道的某个收缩区的空气湍流(类似于噪声)，爆破音是位于声道闭合点处建立的气压及突然地释放。
当激励频率等于震动物体固有的频率时，便以最大振幅来震荡，在该频率上，传递函数有极大值，这种现象称为共振，一个共震体可能存在多个相应强度不同的共振频率。声道是分布参数系统，可以看做是谐振腔，有很多谐振频率。谐振频率由每一瞬间的声道外形决定。这些谐振频率称为共振峰频率，简称共振峰，是声道的重要声学特性。这个线性系统的特征频率特性称为共振峰特性，决定了信号的频谱的总轮廓即包络。为了得到高质量的语音或准确的描述语音，须采用尽可能多的共振峰。在实际应用中，声学语音学中通常考虑前两个峰，语音合成考虑五个共振峰是最现实的。
汉语的特点是音素少，音节少，大约有64个音素，但只有400个左右的音节，即400个左右的基本发音，假如要考虑每个音节有5种音调，也不过有1200多个有调音节即不同的发音。元音属于浊音，脉冲间隔为基音周期，用g(t)表示。其作用于声道，得到的语音信号是g(t)与声道冲激响应h(t)的卷积。g(t)的频谱是间隔为基频的脉冲序列的频谱与声门波频谱[32]的乘积。语音信号可看做便利性随机过程，其统计特性可用信号幅度的概率密度及一些统计量(主要为均值和自相关函数)来描述。对语音的研究表明，其幅度分布有两种近似的形式，较好的为修正Gamma分布[33][34]:

2.2.语音信号的数学模型

表示采样语言信号的离散模型是特别重要的[36]，建立模型的目的是要寻求一种可以表达一定物理状态下的数学关系，而且要使这种关系不仅具有巨大的精度，还要最简单。人们希望模型既是线性的又是时不变的，这是最理想的模型，但根据语音的产生机理，语音信号是一连串的时变过程，不能满足这两种性质。因此我们需要做出一些合理的假设，使得在较短的时间间隔内表示语音信号时，可采用线性时不变模型。在一般的语音信号经典模型中，语音信号被看做线性时不变系统(声道)在随机噪声或准周期脉冲序列下的输出。这一模型用数字滤波器原理进行公式化以后，将称为语音处理技术的基础[37]。
研究表明语音的产生就是声道中的激励，语音传播就是声波在声道中的传播。假若采用流体力学等建立复杂方程的方法进行研究十分复杂。为了简化，通常对声道形状和发音系统进行某些假设，如假设声道是时变的且有不均匀截面的声管，空气流动或声管壁不存在热传导或粘滞消耗，波长大于声道尺寸的声波是沿声管管轴传播的平面波；更进一步简化，进一步假设声道是由半径不同的无损声管级联得到的。在上述这些假设下，得到级联无损声管模型的传输函数，可以证明对大多数语音，该传输函数为全几点函数，只是对鼻音和摩擦音需加入一些零点。但由于任何零点可用多极点逼近，因此可用全极点模型模拟声道[38]。另一方面，级联无损声管与全极点数字滤波器有很多相同的性质，因而用数字滤波器模拟声道特性是一种常用的方法。
产生语音信号的示意图如图2-1所示：

其中N1为斜三角波上升部分的时间，N2为斜三角波下降部分的时间。可见他是一个低通滤波器，通常习惯把它表示为Z变换的全极点模式[41]的形式：

所以整个激励模型可以表示为：

在声管模型中，每个管子可看做一个四端网络，其具有反射系数，这些系数与LPC参数间有唯一的对应关系。声道可由一组截面积或一组反射系数表示。
（2）共振峰模型[43]：
将声道视为谐振腔时，共振峰即为腔体的共振频率。研究表明，用前三个共振峰代表一个元音就可以，而对较复杂的辅音或鼻音，需要用五个以上的共振峰。基于共振峰理论，有三种实用的模型：级联型、并联型和混合型。
级联型认为声道为一组串联的二阶谐振器。根据共振峰理论，整个声道有多个谐振频率和多个反谐振频率(对应声道频率特性的零点)，因而可以被模拟为零极点模型，但对一般元音可用全极点模型，将声道看做一个变截面声管，根据流体力学可得在大多数情况下其为全极点函数，此时共振峰用自回归(AR)模型近似。由于采用LPC技术可以高效的求解AR模型系数[44]，因此该模型应用十分普遍。
对于比较复杂的元音和大部分的辅音，需要采用零极点模型，可用并联型模型表示。但在实际应用中，上述两种模型都较为简单，可用于描述一般的元音，但当鼻化元音或鼻腔参与共振等情况级联模型就不适用了，此时腔体有反谐振特性，需要加入零点，称为极零点模型，此时称为并联型结构。将级联模型和并联模型结合的混合型是较为完备的共振峰模型，其可根据不同性质的语音进行切换。如下图2-3所示：

语音信号模型中，如不考虑周期冲击脉冲串模型E(z)，则斜三角波模型为二阶低通，辐射模型为一阶高通，因而实际信号分析中常采用预加重技术，即对信号取样喉插入一阶高通滤波器，从而只剩下声道部分，便于对声道参数进行分析。

2.3.语音信号的时域、频域特性和短时分析技术

语音信号是随时间变化的非平稳随机过程，因此对于语音的分析一般都是短时分析[46]。这是因为语音虽然是时变的但是具有短时相关性，这个相关性来源于人的发生器官具有惯性，因此语音的状态是不会发生突变，语音在短时间内语音信号的特性基本不变，称之为语音的短时平稳性。语音信号本身是时域信号，对其进行分析时，最直观的方法就是观察其时域波形。语音的短时分析一般需要分帧加窗来处理。语音包含能量这是不难理解的，比如清音和浊音包含的能量就明显不同。
语音信号有许多特点，它的带宽约为5KHz，主要集中在低频段。而且它是一种典型的随机信号：首先，人的每次发声过程都是一个随机过程，很难得到两次完全相同的发音样本；其次，在信号处理中通常假设语音信号是短时平稳的，例如可以认为在语音的浊音段部分，语音的二阶矩统计量是平稳的，即二阶矩平稳，或称为宽平稳。

2.3.1语音信号分帧

为了方便分析输入的音频数据，通常需要进行分帧处理，在分帧中，由于语音信号是是时变的，在短时范围内特征变化较小，所以作为稳态处理，但是超出这个短时范围语音信号就会有变化，因此需要设置相邻两帧间有一部分重叠。在相邻两帧之间基音发生了变化，如正好是两个音节之间或正好是声母向韵母过渡等情况，这时其特征参数有可能变化较大，但是为了使特征参数平滑地变化，在两个不重叠的帧之间插一些帧来提取特征参数，这就形成了相邻帧之间有重叠部分[47]。
设输入的音频文件长度为N，取每帧长为wlen。为了使其平滑过渡，两帧之间需要插入一帧或几帧，后一帧相对前一帧的位移量用inc表示，则相邻两帧之间的重叠部分为：

因此，每一帧在音频中的位置为：

短时能量与短时平均幅度函数的主要用途有：区分浊音段与清音段，因为浊音的短时能量比清音大得多；区分声母与韵母的分界；区分无话段与有话段的分界。
（2）短时平均过零率：
短时平均过零率是语音信号时域分析中的一种特征参数。它是指每帧内信号通过零值的次数。对有时间横轴的连续语音信号，可以观察到语音的时域波形通过横轴的情况。在离散时间语音信号情况下，如果相邻的采样具有不同的代数符号就称为发生了过零，因此可以计算过零的次数。单位时间内过零的次数就称为过零率。一段长时间内的过零率称为平均过零率。如果是正弦信号，其平均过零率就是信号频率的两倍除以采样频率，而采样频率是固定的。因此过零率在一定程度上可以反映信号的频率信息。
设语音波形时域信号为x(n)分帧处理后得到第i帧语音信号为yi(n):，帧长为L，短时平均过零率为：

选择的窗函数和窗宽的不同，对短时傅里叶谱的影响是不同的。
可以看出在矩形窗和汉明窗两种窗函数下，短时频谱图都有两种变化：由周期性激励引起的快变化，反映了基音频率的各次谐波；由声道的共振峰引起的慢变化，反映了各共振峰的频率和带宽。还可以看出两个频谱图之间存在明显的差别。采用矩形窗时，基音谐波的各个峰都比较尖锐，且整个频谱图显得比较破碎（类似于噪声），这是因为矩形窗的主瓣较窄，具有较高的频率分辨率，但它也具有较高的旁瓣，因而使基音的相邻谐波之间的干扰比较严重。在相邻谐波间隔内有时叠加，有时抵销，出现了一种随机变化的现象。相邻谐波之间的这种严重“泄露”现象，抵消了矩形窗主瓣窄的优点，因此，在语音短时频谱分析中极少采用矩形窗。当加汉明窗时，得到的短时频谱要平滑的多，因而在语音分析中汉明窗用得比较普遍。
综上所述，关于短时谱和移动窗可以得到以下结论。
1.长窗具有较高的频率分辨率，较低的时间分辨率。从一个基音周期到另一个基音周期，共振峰是要发生变化的，这一点即使从语音波形上也能够看出来。然而如果采用较长的窗，这种变化便被模糊了，因为长窗起到了时间上的平均作用。
2.短窗具有较低的频率分辨率，较高的时间分辨率。采用矩形窗时，能够从短时谱中提取出共振峰从一个基音周期到另一个基音周期所发生的变化。当然，激励源的谐波的细致结构也从短时频谱图上消失了。
3.窗宽的选择需折中考虑。短窗具有较好的时间分辨率，能够提取出语音信号中的短时变化，损失了频率分辨率。但应注意到，语音信号的基音周期提取范围很大。因此，窗宽的选择应当考虑到这个因素。
4.矩形窗和汉明窗的频谱特性都具有低通的性质，在截止频率出都比较尖锐，当其通带都比较窄时（窗越宽，其带通越窄），加窗后得到的频谱能够很好逼近短时语音信号的频谱。窗越宽，逼近效果越好。

2.4.带噪语音处理

2.4.1带噪语音的产生

噪音在生活中无处不在，我们人耳其实已经习惯于噪音的环境，日常生活里绝对没有噪音的环境不存在。噪声主要来于两个方面：一方面是设备噪声，声音进入麦克风是模拟信号，通过声卡转为数字型号，数字信号是计算机可以处理的信号，处理结束又转换为模拟信号输出，也就是耳机音箱最终播放出的声音。也就是说录一段音频其实要经过好多设备处理，在这个过程中音频难免就会受到设备电路的干扰；另一方面是地噪，严格讲地噪也是设备噪音，但是不是设备自身造成的，而是设备外部电路造成，比如家庭电路没有接地，电流干扰到了麦克风、声卡等电路，所以有些声卡带有专门的接地端，就是为了避免地噪的干扰。
（1）信噪比：
对于带噪信号而言，其所带的噪声“量”如何衡量是极其重要的，这时我们就需要通过信噪比来表示噪声的能量信噪比是信号和噪声能比值的对数值。它被定义为：

高斯白噪声的功率谱密度服从均匀分布，幅度分布服从高斯分布，在任意两个不同时刻上的随机变量之间，不仅是互不相关的，而且还是统计独立的。

2.4.2小波变换去噪

小波去噪是通过短波实现噪音消除[52]，与高斯去噪的基本原理一致。近年来，小波理论得到了非常迅速的发展，而且由于其具备良好的时频特性，因而实际应用也非常广泛。在去噪领域中，小波理论也同样受到了许多学者的重视，他们应用小波进行去噪并获得了非常好的效果。
小波去噪方法包括三个基本的步骤：对含噪声信号进行小波变换；对变换得到的小波系数进行某种处理，以去除其中包含的噪声；对处理后的小波系数进行小波逆变换，得到去噪后的信号。小波去噪方法的不同之处集中在第一步。小波去噪实际上是特征提取和低通滤波功能的综合，其流程框图如图2-4所示：

【下面这两页粘贴后内容缺失较多，直接放截图了】

3 matlab程序

（非高斯白噪声，与说明书内容略有不符，读者可查阅其他博客做简单更改）

function varargout = yuyinshuzhixinhao(varargin)% YUYINSHUZHIXINHAO MATLAB code for yuyinshuzhixinhao.fig%      YUYINSHUZHIXINHAO, by itself, creates a new YUYINSHUZHIXINHAO or raises the existing%      singleton*.%%      H = YUYINSHUZHIXINHAO returns the handle to a new YUYINSHUZHIXINHAO or the handle to%      the existing singleton*.%%      YUYINSHUZHIXINHAO('CALLBACK',hObject,eventData,handles,...) calls the local%      function named CALLBACK in YUYINSHUZHIXINHAO.M with the given input arguments.%%      YUYINSHUZHIXINHAO('Property','Value',...) creates a new YUYINSHUZHIXINHAO or raises the%      existing singleton*.  Starting from the left, property value pairs are%      applied to the GUI before yuyinshuzhixinhao_OpeningFcn gets called.  An%      unrecognized property name or invalid value makes property application%      stop.  All inputs are passed to yuyinshuzhixinhao_OpeningFcn via varargin.%%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one%      instance to run (singleton)".%% See also: GUIDE, GUIDATA, GUIHANDLES% Edit the above text to modify the response to help yuyinshuzhixinhao% Last Modified by GUIDE v2.5 18-Jun-2021 14:58:02% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name',       mfilename, ...    'gui_Singleton',  gui_Singleton, ...    'gui_OpeningFcn', @yuyinshuzhixinhao_OpeningFcn, ...    'gui_OutputFcn',  @yuyinshuzhixinhao_OutputFcn, ...    'gui_LayoutFcn',  [] , ...    'gui_Callback',   []);if nargin && ischar(varargin{1})    gui_State.gui_Callback = str2func(varargin{1});endif nargout    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});else    gui_mainfcn(gui_State, varargin{:});end% End initialization code - DO NOT EDIT% --- Executes just before yuyinshuzhixinhao is made visible.function yuyinshuzhixinhao_OpeningFcn(hObject, eventdata, handles, varargin)% This function has no output args, see OutputFcn.% hObject    handle to figure% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)% varargin   command line arguments to yuyinshuzhixinhao (see VARARGIN)% Choose default command line output for yuyinshuzhixinhaohandles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes yuyinshuzhixinhao wait for user response (see UIRESUME)% uiwait(handles.figure1);% --- Outputs from this function are returned to the command line.function varargout = yuyinshuzhixinhao_OutputFcn(hObject, eventdata, handles) % varargout  cell array for returning output args (see VARARGOUT);% hObject    handle to figure% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)% Get default command line output from handles structurevarargout{1} = handles.output;% ##############################################################%           录音%###############################################################% --- Executes on button press in pushbutton1.function pushbutton1_Callback(hObject, eventdata, handles)% hObject    handle to pushbutton1 (see GCBO)% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)recObj = audiorecorder;disp('Start speaking.')%recordblocking(recObj, 5);recordblocking(recObj,5);disp('End of Recording.');% 回放录音数据axes(handles.axes1);%表示的是将上面的坐标轴做为当前坐标轴,在其上做图.play(recObj);% 获取录音数据myRecording = getaudiodata(recObj);% 绘制录音数据波形plot(myRecording);%存储语音信号filename = './data1.wav'; audiowrite(filename,myRecording,8000);% ##############################################################%           FFT频谱图%###############################################################% --- Executes on button press in pushbutton2.function pushbutton2_Callback(hObject, eventdata, handles)% hObject    handle to pushbutton2 (see GCBO)% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)fs=22050;% 读取音频文件filename=strcat('data1.wav');[x,fs]=audioread(filename);ainfo = audioinfo(filename);bits = ainfo.BitsPerSample;p=audioplayer(x,fs,bits);sound(x,fs,bits);y1=fft(x,1024);%对信号做1024点FFT变换f=fs*(0:511)/1024;%--------------------------[H,w]=freqz(x) ;%绘制原始语音信号的频率响应图Hf=abs(H);  %取幅度值实部Hx=angle(H);  %取相位值对应相位角%clf% %-----------------------------------global h1 h2 h3 h4;h=0;if ishandle(h1)  %判断h1是否为句柄    delete(h1);    h=1;endif ishandle(h2)  %判断h2是否为句柄    delete(h2);    h=1;endif ishandle(h3)  %判断h3是否为句柄    delete(h3);    h=1;endif ishandle(h4)  %判断h4是否为句柄    delete(h4);    h=1;endif h    axes('parent',handles.uipanel2) %重建一个axesend%------------------------------------------------%axes(handles.axes1);%axes('parent',handles.uipanel3) subplot(221);plot(x);title('原始语音信号时域图') xlabel('时间/S');ylabel('幅度');%figure(1)%axes(handles.axes1);%axes('parent',handles.uipanel2)  %重建一个axes%绘制频率响应图%--------------------------------------------% %subplot(2,1,1);% % axes(handles.axes1);% % axes('parent',handles.uipanel1) subplot(222);plot(w,20*log(Hf)); %幅值变换为分贝单位title('离散系统幅频特性曲线') xlabel('频率/HZ');ylabel('幅度');% % %----------------------% axes(handles.axes4);% axes('parent',handles.uipanel6) subplot(223);plot(w,Hx)title('离散系统相频特性曲线')xlabel('频率/HZ');ylabel('相角');% %figure(1) % %freqz(x);% %title('频率响应图');% %---------------------------------------------% %原始信号频谱% axes(handles.axes3);% axes('parent',handles.uipanel5) % %figure(2)subplot(224);plot(f,abs(y1(1:512)));%title('原始信号频谱','color','b');title('原始语音信号频谱');xlabel('频率/HZ');ylabel('幅度/dB');% ##############################################################%           叠加噪声%###############################################################% --- Executes on button press in pushbutton3.function pushbutton3_Callback(hObject, eventdata, handles)% hObject    handle to pushbutton3 (see GCBO)% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)%fs=22050;% 读取音频文件filename=strcat('data1.wav');[x,fs]=audioread(filename);ainfo = audioinfo(filename);bits = ainfo.BitsPerSample;p=audioplayer(x,fs,bits);sound(x,fs,bits);y1=fft(x,1024); %对信号做1024点FFT变换f=fs*(0:511)/1024;% %-------------------------------------------------------------%加噪声L=length(x);      %计算音频信号的长度%x1=rand(1,length(x)); %产生与x长度一致的随机信号x1=0.1*randn(L,2);x2=x1+x;sound(x2);y2=fft(x2,1024);% %------------------------------------------------------------global h1 h2 h3 h4;h=0;if ishandle(h1)  %判断h1是否为句柄    delete(h1);    h=1;endif ishandle(h2)  %判断h2是否为句柄    delete(h2);    h=1;endif ishandle(h3)  %判断h3是否为句柄    delete(h3);    h=1;endif ishandle(h4)  %判断h4是否为句柄    delete(h4);    h=1;endif h    axes('parent',handles.uipanel2) %重建一个axesend%------------------------------------------------%axes(handles.axes1);subplot(221);plot(x);title('原始语音信号时域图') % %-----------------------------subplot(222);plot(x2);title('加高斯噪声后的语音信号时域图') xlabel('时间');ylabel('幅度');%-----------------------------subplot(223);plot(abs(y1));title('原始语音信号频谱') xlabel('HZ');ylabel('幅度');% %-----------------------------subplot(224);plot(abs(y2));title('加噪声后的语音信号频谱') xlabel('HZ');ylabel('幅度');% ##############################################################%           butterworth模拟滤波器%###############################################################% --- Executes on button press in pushbutton4.function pushbutton4_Callback(hObject, eventdata, handles)% hObject    handle to pushbutton4 (see GCBO)% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)% 双线变换法设计巴特沃斯低通滤波器% 读取音频文件filename=strcat('data1.wav');[x,fs]=audioread(filename);ainfo = audioinfo(filename);bits = ainfo.BitsPerSample;p=audioplayer(x,fs,bits);sound(x,fs,bits);y1=fft(x,1024); %对信号做1024点FFT变换f=fs*(0:511)/1024;% %-------------------------------------------------------------%加噪声L=length(x);      %计算音频信号的长度%x1=rand(1,length(x)); %产生与x长度一致的随机信号x1=0.1*randn(L,2);x2=x1+x;sound(x2);y2=fft(x2,1024);  %滤波处理前%--------------------------------------------------------------wp=0.1*pi;disp(wp)ws=0.4*pi;disp(ws)Rp=1;Rs=15;Fs=22050;Ts=1/Fs;%将模拟指标转换成数字指标wp1=2/Ts*tan(wp/2);  ws1=2/Ts*tan(ws/2);%选择滤波器的最小阶数%[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %[N,Wn]=cheb2ord(wp1,ws1,Rp,Rs,'s'); %切比雪夫2[N,Wn]=ellipord(wp1,ws1,Rp,Rs,'s'); %椭圆disp(N)disp(Wn)%创建butter worth模拟滤波器[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);%用双线形变换法实现模拟滤波器到数字滤波器的转换[bz,az]=bilinear(b,a,Fs);f1=filter(bz,az,x2);%----------------------------------------------------%绘制频率响应曲线[H,W]=freqz(bz,az);%grid%----------------------------------------------------global h1 h2 h3 h4;h=0;if ishandle(h1)  %判断h1是否为句柄    delete(h1);    h=1;endif ishandle(h2)  %判断h2是否为句柄    delete(h2);    h=1;endif ishandle(h3)  %判断h3是否为句柄    delete(h3);    h=1;endif ishandle(h4)  %判断h4是否为句柄    delete(h4);    h=1;endif h    axes('parent',handles.uipanel2) %重建一个axesend%------------------------------------------------%axes(handles.axes1);subplot(2,2,1);plot(x2);title('滤波前的时域波形')%-------------------------subplot(2,2,2);plot(f1);title('滤波后的时域波形')%-------------------------%播放滤波后的信号sound(f1);F0=fft(f1,1024); %对信号做1024点FFT变换f2=fs*(0:511)/1024;%-------------------------subplot(2,2,3);plot(f2,abs(y2(1:512)));title('滤波前的频谱')xlabel('HZ');ylabel('幅度');%-------------------------subplot(2,2,4);plot(f2,abs(F0(1:512)));title('滤波后的频谱')xlabel('HZ');y
来源：Unite One
                                                        声明：本站部分文章及图片转载于互联网，内容版权归原作者所有，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！