DSP期末总复习
 
基于教材：《数字信号处理——原理、实现及应用（第3版）》 高西全等 编著
 
一、基本概念
 
1.1 信号
 
1.2 系统
 
线性性
时不变性
因果性/可实现性
 （1）def: 系统n时刻的输入取决于n时刻和n时刻以前的输入信号，和n时刻之后的输入信号无关
 （2）充要条件：n<0时h(n)=0
 （3）因果系统的单位脉冲响应必然是因果序列
稳定性
 充要条件：系统的单位脉冲响应绝对可和
 
1.4 差分方程
 
def：从0开始，到某项结束
递推法求解
 
二、傅里叶变换
 
2.1 连续信号：CTFT和CTFS
 
记住公式，不做证明。
 
 
2.2 时域离散信号：DTFT和DTFS、周期信号的傅里叶变换
 
2.2.1 DTFS
 
公式、证明如下：
 
频谱是离散线状谱，周期为N.（虽然横坐标不是频率，但是是对频率的采样，所以也是频谱）
 【变量在n,k之间切换】这实质上体现了：离散谱是对包络线的采样。
 
2.2.2 DTFT
 
适用于非周期信号，是通过令周期序列的周期趋于无穷得到的，同时还要求序列绝对可和。
 公式：
 
【变量在n,w之间切换】频谱是连续的，是对DTFS离散谱的包络，周期为2pi.
 
2.2.3 周期信号的傅里叶变换
 
公式及证明：
 
 
2.3 傅里叶变换的性质和常用结论
 
性质
 （1）周期性
 （2）线性性
 （3）移位性：x(n-n0) --> 外面乘一个e^-jwn0
 （4）频移性：e^-jw0n --> X里面变成w-w0
 （5）频域微分性质
 （6）对称性
 （7）卷积定理
基本序列的FT
 
 
2.4 离散傅里叶变换（DFT）
 
DTFT和ZT的变换结果都是连续函数，而DFT既是时域到频域的变换，其变换结果又是对时域离散信号频谱的等间隔采样。
 
2.4.1 基本概念
 
设序列长度为M，定义N点DFT为：
 
 要求：N>=M，注意WN的代换
 
2.4.2 关系
 
N点DFT是序列的Z变换在单位圆上的等间隔采样
N点DFT是序列的DTFT在[0,2pi]上的等间隔采样
IDFT是DTFS的主值区序列，X(k)是X(k)~的主值序列
 
2.4.3 性质
 
线性性
隐含周期性
有限长序列的循环移位：先周期延拓，再左移m个单位，最后取主值序列。
 循环移位性质（证明见下）
复共轭序列的DFT
实信号DFT的特点
 
 
2.5 卷积
 
2.5.1 线性卷积
 
定义
 运算过程：将x(n)中的n换成m，h(n)中的n换成n-m，相同m的序列值相乘后，再对乘积求和。
图解法
解析法【注意分类讨论】
性质：移位、运算律
 
2.5.2 循环卷积
 
三个长度量：L：循环卷积的长度 M：x(n)的长度 N：h(n)的长度
 x((n-m))L是以L为周期的信号。
 形成循环倒相序列，为了解决掉负数项，由周期性给系数加L，n值由0取到L-1，得到矩阵
 
2.6 快速傅里叶变换（FFT）
 
2.6.1 如何减小运算量？
 
不断将长序列的DFT分解为短序列的DFT，并利用周期性和对称性及一些特殊值来减少DFT运算量。
 
2.6.2 基2 FFT算法
 
按n的奇偶性分解
引入l,进行变量代换
利用WN的性质，两个角标同时除以2，得到两个N/2点DFT
得到周期性公式
蝶形图：（以8点为例）
 （1）最后求得的八个数，是两两一组的形式求出来的，即X(k),X(k+N/2)
 为了解出这两个未知数，我们需要知道对应的X1(k),X2(k)
 （2）X1(k)又和X1(k+N/4)构成另一对，可以看出，两组的间距、运算点数每次除以2
 ……
 （3）直到到达最后一层
 WN的角标变化：两组[0 1 2 3] N --> 四组[0 1] N/2【实际写成[0 2] N】 --> 八组[0] N/4【实际写成[0] N】
 
2.7 模拟信号的频谱分析
 
由模拟信号时域采样得到N点采样序列，经过DFT，再乘以T，就是模拟信号在频域的采样。
 f和k可以进行坐标转换，以f为横坐标，以X(k)为函数值，绘制出|X(k)|的包络/X(k)的相位包络，就是模拟信号的幅频/相频特性。
参数选择
 
误差分析
 连续信号的谱分析主要关心两个指标：分辨率、谱分析范围
 （1）频谱混叠：若采样频率不满足采样定理，会在0.5Fs附近引起混叠
 （2）截断效应：CTFT是在无穷区间上的，然而实际观察到的模拟信号一般是有限长的，没有观察到的是0.
 用矩形窗截断。
 1）频谱泄漏：原来离散的谱线向两边展宽。矩形窗长度越长，展宽的宽度越窄。
 泄漏会使频谱模糊，谱的分辨率降低。
 2）谱间干扰：矩形窗函数的频谱存在许多旁瓣，和主信号频谱卷积后形成很多旁瓣，旁瓣起着谱间干扰的作用。
 （3）栅栏效应：N点DFT得到的只是N个采样点上的频谱值，两点间的频谱是不知道的，就像被栅栏遮住一样。
DFT对周期信号进行谱分析
 
三、Z变换
 
3.1 Z变换
 
def: 不加说明均指双边Z变换
 只有幂级数收敛，Z变换才有意义–>收敛域
Z变换和傅里叶变换的关系：傅里叶变换是单位圆上的Z变换（前提是Z变换必须存在）
四种序列的收敛域
 因果序列的收敛域：R<|z|<=INF 零极相消不影响收敛！
常用结论证明：
 
 
 
3.2 Z变换的性质
 
线性性：注意收敛域取交集！
序列移位：x(n-n0) --> 提取出z^(-n0)
时间反转：n取负 --> z取倒数
z域微分：nx(n) --> -z*(X(z))’
共轭序列：对x(n)取共轭 --> X(z)同时对z和X取共轭
五个定理：
 （1）时域卷积定理
 
3.3 逆Z变换
 
部分分式展开法
 最关键的就是最后一步时，根据收敛域判别可能的原序列！
 以及，是否足够熟练地作出准确记忆和判断。【助记：z大于啥的时候正常，小于加负号】
围线积分法
 其中：c是X(z)的收敛域中一条包含原点的逆时针旋转的封闭曲线，zk是F(z)在围线c内的M个极点
 
 进行留数相关运算时，z都写成z-a的形式，不要带着z^-1的形式了。
例题
 （1）课本2.3.6【要点：F(z)的极点与n的取值有关；因果序列的性质；多阶极点的处理】
 
 （2）课本2.3.10【要点：z变换性质的综合应用】
 
 （3）课后2-24【要点：按收敛域分类讨论】
 
 
3.4 Z变换的应用
 
传输函数/频率响应函数、系统函数
根据极点分布，分析系统的因果性、稳定性
 （1）因果系统的系统函数的极点均在某个圆内，收敛域包含无穷
 （2）系统稳定时，系统函数收敛域一定包含单位圆/系统函数的极点不能位于单位圆上。
 （3）因果稳定的条件：H(z)的极点集中在单位圆内
根据零极点分布，分析系统的频率特性
 （1）基本概念
 设：M个零点，N个极点
 频率w由0增大到2pi的过程中，OB矢量绕坐标原点逆时针转一圈。
 极点零点分别指向B --> 极点矢量，零点矢量
 幅频特性 = 零点矢量长度之积/极点矢量长度之积
 相频特性 = w(N - M) + 零点矢量相角之和 - 极点矢量相角之和
 （2）零极点分布对幅频特性的影响：
 系统函数的极点主要影响幅频特性的峰值，峰值频率在极点附近。极点越靠近单位圆，峰值越高越尖锐。若极点在单位圆上，幅度无穷大，系统不稳定。
 系统函数的零点主要影响幅频特性的谷值，谷值频率在零点附近。零点越靠近单位圆，谷值越接近零。若零点在单位圆上，谷值为零。
 处于坐标原点的零极点不影响幅频特性。
 （3）例题：课本2.4.6
 
 （4）例题：课本2.4.8
 
 
四、模拟信号数字处理
 
这一章要干什么？
 
用数字信号处理方法处理模拟信号。处理思路如下：
 预滤波 --> ADC --> 数字信号处理 --> DAC --> 平滑滤波
 
4.1 预滤波
 
采样前对模拟信号进行预滤波，预滤波器就是一个模拟低通滤波器，阻带截止频率0.5 \omega s
 也称抗混叠滤波器。
 
4.2 采样
 
4.2.1 理想采样信号
 
时域【横坐标都是t】：模拟信号 + 单位冲激串调制 —> 理想采样信号
频域【横坐标都是模拟角频率】【CTFT得到】：理想采样信号的频谱是模拟信号频谱以采样频率为周期进行周期延拓。
原理推导：
 
 
4.2.2 时域离散信号
 
区别：模拟信号和理想采样信号在整个t轴上有定义，而时域离散信号仅在整数点有定义。
时域【横坐标是n】：模拟信号 + 等间隔T采样 —> 时域离散信号
频域【横坐标是数字频率】【DTFT得到】：时域离散信号的频谱是模拟信号频谱变换一下横坐标（除以T）。
采样频率必须大于等于模拟信号最高频率的两倍以上，才能不丢失信息。
原理推导：
 
 
4.4 理想恢复
 
def: 如何用理想采样信号恢复模拟信号？
原理推导：将理想采样信号通过理想低通滤波器G
 
 g(t)的作用：在采样点之间连续插值
 
4.5 线性模拟系统的数字模拟
 
用数字系统模拟线性模拟系统的外部特性
传输函数的关系推导：
 
 
五、滤波器
 
5.1 一些散落的实例
 
5.1.1 滑动平均滤波器
 
5.1.2 全通滤波器
 
5.1.3 最小相位滤波器
 
5.2 模拟滤波器设计
 
各种滤波器都是基于低通滤波器设计。
 
5.2.1 设计指标
 
模拟滤波器的设计流程：
 （1）根据信号处理要求确定设计指标
 （2）选择滤波器类型
 （3）计算滤波器设计
 （4）确定滤波器系统函数
 （5）综合实现
因为几种典型的AF的相位特性已知，所以一般只给出对幅频特性的要求，而不考虑相位特性（若有要求则需另外考虑）。
 系统函数：Ha(s) 频率响应函数/传输函数：Ha(j/omega)
 通带边界频率：/omega p 阻带边界频率：/omega s
 通带最大衰减：/alpha p 阻带最小衰减：/alpha s 【注意是dB】
 损耗/衰减函数 3dB截止频率
 
 
5.2.2 巴特沃斯滤波器
 
 
5.2.3 切比雪夫滤波器
 
I型滤波器：幅频特性在通带等波纹，阻带单调下降
 II型滤波器：幅频特性在阻带等波纹，通带单调下降
 
 
5.2.4 椭圆滤波器
 
幅频特性在通带和阻带都等波纹。通带和阻带波纹固定时，阶数越高，过渡带越窄；阶数固定时，通带和阻带波纹越小，过渡带越宽。
 阶数相同时，椭圆滤波器可以获得对理想滤波器幅频响应的最好逼近。
 
5.3 数字滤波器设计
 
5.3.1 基本概念
 
数字滤波：通过对输入信号进行数值运算，让输入信号中有用的频率成分以较高的保真度通过，滤除某些无用的频率成分，实现对输入信号的选频处理。本章主要介绍选频性滤波器。
数字滤波器的优点：精度高，稳定性好，体积小，实现方法灵活，不存在阻抗匹配问题
分类：无限长单位脉冲响应数字滤波器IIRDF，有限长单位脉冲响应数字滤波器FIRDF
 系统函数分别为：
 
设计方法：间接设计法【先根据数字滤波器指标设计相应的过渡模拟滤波器，再转换成DF】、直接设计法
理想滤波器的特点
 （1）通带内幅度为常数
 （2）阻带中幅度为0
 （3）单位脉冲响应是非因果无限长序列
 （4）物理不可实现
群时延：相频特性的负导数 群时延为常数表示输入信号中的所有频率分量的时间延迟相同。
信号通过理想滤波器波形不失真的充要条件：通带内幅度特性是常数，并且具有线性相位
理想滤波器的近似实现
 （1）h(n)向右移动，忽略n<0的较小的部分
 （2）h(n)对中心最大值对称则具有线性相位特性
设计实质：根据设计指标确定一个因果稳定的系统函数H(z)
 
5.3.2 IIRDF
 
一、基本概念
 
间接法设计流程：
 （1）给定数字滤波器指标
 （2）将数字滤波器指标转换成相应的过渡模拟滤波器指标
 （3）选择合适类型的模拟滤波器
 （4）设计满足要求的过渡模拟滤波器Ha(s)
 （5）将Ha(s)转换成数字滤波器H(z)
转换：幅频特性函数、三个带、频率
 
 转换实质：用一种从s平面到z平面的映射函数将Ha(s)转换成H(z)
 映射函数的要求：
 （1）因果稳定的模拟滤波器转换后仍因果稳定
 （2）H(z)的频率响应特性能近似模仿Ha(s)的片断常数频率响应特性
 
二、脉冲响应不变法（这两种方法都是转换方法）
 
基本思想：【波形逼近】（设模拟滤波器的单位冲激响应为ha(t)）对ha(t)等间隔采样，得到数字滤波器的单位脉冲响应h(n)
原理推导、性能分析
 
 模拟系统因果稳定，Ha(s)的所有极点位于s平面的左半平面，这些极点全部映射到z平面单位圆内，因此数字滤波器也因果稳定。
 数字滤波器频率响应是模拟滤波器频率响应的周期延拓。
优点：（若模拟滤波器具有带限特性【实际一般都不满足】，T满足采样定理）数字滤波器频率响应完全模仿了模拟滤波器频率响应。
 缺点：总是存在频谱混叠失真，混叠严重时不满足阻带衰减指标，不适合设计高通和带阻滤波器。
 如何减小频谱混叠失真？
 （1）选用具有锐截止特性【锐截止特性：滤波器从通带能够很快的过渡到阻带的滤波器特性，即过渡带很窄】的模拟滤波器【代价：阶数升高】
 （2）提高采样频率【代价：对处理速度要求更高】
例题6.2.1【一般取T = 1s】
 
三、双线性变换法
 
基本思想：【算法逼近】一般数字滤波器用差分方程描述，模拟滤波器用微分方程描述。因此只要能用差分近似微分，就可以将差分方程转换成微分方程。
原理推导、性能分析
 
 s平面的整个正虚轴映射成有限宽度的数字频段[0,pi]，引起了数字频率和模拟频率之间的严重非线性畸变。
优点：从原理上彻底消除了频谱混叠
 缺点：仅能保持原模拟滤波器的片断常数幅频响应特性，不能保持相频响应特性。
设计步骤
 （1）确定数字滤波器四个指标。
 （2）非线性预畸变校正，转换成相应的模拟滤波器指标。
 （3）设计满足指标要求的过渡模拟滤波器。
 （4）由映射关系，转换成数字滤波器
采样间隔T的选择
 脉冲响应不变法：T = 1s 双线性变换法：T = 2s
例题6.2.3（略）
 
5.3.3 FIRDF
 
最大优点：可以实现线性相位滤波
 其他优点：全零点滤波器，硬件和软件实现结构简单，不用考虑稳定性问题。
 当幅频特性相同时，FIRDF阶数比IIRDF高得多。
 
一、线性相位FIRDF
 
1. 线性相位
 重新表示频率响应函数：幅度特性函数相位特性函数 【区别：幅频响应函数，相频响应函数】
 （1）第一类线性相位FIRDF的相位特性函数：/theta(w) = -w /tao
 （2）第二类线性相位FIRDF的相位特性函数：/theta(w) = /theta0 - w /tao
 第二类线性相位FIRDF只讨论theta0 = -pi/2的情况。
 线性相位滤波器又称恒定群延时滤波器。
 2. 时域约束（对h(n)的约束）
 第一类：h(n)关于（N-1）/2偶对称
 第二类：h(n)关于（N-1）/2奇对称
 3. 频域约束
 （1）h(n) = h(N - 1 - n) N为奇数，Hg(w)关于0 pi 2pi三点偶对称：可以实现低通、高通、带通、带阻
 （2）h(n) = h(N - 1 - n) N为偶数，Hg(w)关于0 2pi两点偶对称，关于pi奇对称 ：可以实现低通、带通
 （3）h(n) = - h(N - 1 - n) N为奇数，Hg(w)关于0 pi 2pi三点奇对称：可以实现带通
 （4）h(n) = - h(N - 1 - n) N为偶数，Hg(w)关于0 2pi两点奇对称，关于pi偶对称：可以实现高通、带通
 4. 零点分布
 若zk是零点，则zk^-1,zk, zk* ^-1都是零点，四个一组
 
二、窗函数法设计
 
**1. 基本思想：**用FIRDF逼近希望的滤波特性
 窗函数设计法：截取hd(n)有限长的一段因果序列，并用合适的窗函数进行加权作为FIRDF的单位脉冲响应h(n)
 2. 基本方法
 （1）确定技术指标
 （2）构造希望逼近的频率响应函数
 （3）求hd(n)
 （4）加窗得到FIRDF的单位脉冲响应h(n)
 
 3. 设计性能分析
 （1）理论推导
 逼近误差实际就是加窗的影响，窗函数的类型、长度都会影响逼近误差。
 
 （2）两点误差（加窗效应）
 1）理想低通滤波器过渡带宽度为0，但是：Hg(w)以wc为中心形成过渡带，过渡带宽度 = 主带宽度
 2）理想低通滤波器通带阻带幅度为1,0，但是：Hg(w)在通带阻带均有波纹。旁瓣越大，波纹幅度越大。
 （3）一般的设计指标会给出：过渡带宽度、通带最大衰减、阻带最小衰减
 要想使过渡带变窄：增大窗函数长度N
 要想满足衰减：改进窗函数形状
 4. 典型窗函数
 （1）矩形窗
 三个参数：旁瓣峰值、过渡带宽度、阻带最小衰减
 （2）汉明窗
 （3）凯塞窗
 1）一种可调整的窗函数，最有用且近似最优。通过调整控制参数可以达到不同的阻带最小衰减，提供最小的主瓣宽度，也就是最窄的过渡带。
 2）参数公式
 5. 基于窗函数的设计流程
 （1）给定技术指标四个
 （2）选择窗函数
 （3）确定窗口参数N,beta
 （4）计算截止频率wc
 （5）确定理想滤波器：离散傅里叶反变换 --> hd[n]
 （6）对称加窗截断，得到滤波器
 
5.3.4 IIR数字滤波器与FIR数字滤器的典型区别
 
性能上：
 ●IIR经济高效，可用较低的阶数获得较好的选择性。【IIR滤波器系统函数的极点可以位于单位圆内的任何地方】但是选择性越好，相位非线性越严重。
 ●FIR滤波器可以得到严格的线性相位特性，但成本较高，只能用较高的阶数达到高的选择性。【FIR滤波器系统函数的极点固定在原点】
结构上：
 ●IIR必须采用递归结构：极点位置必须在单位圆内，否则系统不稳定。
 ●FIR主要采用递归结构：不存在稳定问题。可以采用FFT实现，在相同阶数条件下，运算速度可以大大提高。
设计工具上：
 ●IIR滤波器可以借助模拟滤波器的成果，一般都提供有效的封闭形式的设计公式，计算工作量比较小，对计算工具要求不高。
 ●FIR滤波器一般没有有效的封闭形式的设计公式，对计算工具要求较高。
应用上：
 ●IIR滤波器主要用于设计具有片断常数特性的滤波器【低通、高通、带通、带阻】。
 ●FIR滤波器能适应某些特殊的应用。
 
六、时域离散系统的实现
 
6.1 基本概念
 
基本实现方法：软件实现、硬件实现
本章重点研究差分方程/系统函数的实现方法
一些指标
 
6.2 FIR网络结构
 
FIR的单位脉冲响应有限长，差分方程和系统函数如下：
 
 
直接型结构：长度为N的FIRDF直接型结构需要N-1个单位延时器，N个乘法器，N-1个加法器
级联型结构
 （1）先因式分解，整理成一阶二阶网络的级联形式
 （2）其中每个网络都用直接型结构实现
线性相位结构
 公式：
 
 若N为偶数，乘法器减少到0.5N个；若为奇数，乘法器减少到0.5（N+1）个
 
6.3 IIR网络结构
 
特点：信号流图中含有反馈支路，其单位脉冲响应无限长。
 
 
直接型结构
 
 
级联型结构
 
 
并联型结构
 
 
转置型结构
 
 
附录：利用复数、级数知识证明定理
 
1.DTFT的时域卷积定理
 
 
2.DTFT的频域卷积定理
 
 
3.DTFT的巴塞伐尔定理
 
 
4.ZT的时域卷积定理
 
 
5.ZT的频域卷积定理
 
 
6.ZT的巴塞伐尔定理
 
 
7.DFT的时域循环卷积定理
 
 
8.DFT的频域循环卷积定理
 
 
9.离散巴塞伐尔定理
 

试卷框架
 
一、填空题
 10个空，每个2分，10x2=20
 二、计算
 三、设计
 四、问答
 一题（5分）至多两道
 
复习要点
 

 
 
 
 
 
 
 
复习知识点
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
考点分类
 
1.IIR低通数字滤波器设计
 
脉冲响应不变法
 
 
 
 
 
 
 
 
 
 
 
 
 双线性变换步骤
 
 
 
 
 
1.巴特沃思低通滤波器设计思路，基本步骤
 1.1确定数字性能指标，频率指标
 1.2求出模拟域的性能指标，边缘频率指标
 1.3根据巴特沃思型幅度平方函数求模拟滤波器阶数N
 1.4求通带截止频率wc
 1.5求归一化巴型滤波器模拟域表征系统函数Han(s)
 1.6由双线性变换把Han(s)变换到数字域系统函数h(z)
 
2.双线性设计原理及优缺点
 原理 数字滤波器频率响应和模拟滤波器频率响应做相似变换，使得角频率之间是单极映射，可以避免频率响应的混叠失真，其中omg=2/T tan(w/2)，s=2/T*(z-1)/(z+1)
 优点 可以避免频率响应的混叠失真
 缺点 频率之间的变换是非线性的，计算复杂度高
 A
 B
 C
 步骤要清晰
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
 
 
2.FIR低通数字滤波器设计
 

 
 
 
 
 
 
 
 
 各种窗函数
 
 
 
 
 
 
 设计
 
 
 
 
①理想Hd(ejw)→hd(n)
 ②截短hd(n)
 h(n)=hd(n)w(n) 窗函数
 N 合适窗函数
 ③h(n)→H(ejw)→Hd(ejw)
 
例题
 

 
 
 
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
3.FFT 8点 DIT 原理 信号流图 （8分）
 

 
 

 
 
 
 
 
 
 
 
 
 
 
 
4.DFT 与 DTFT 关系及z变换关系
 

 
 几大公式 物理含义
 w → k
 DFT和DTFT关系：
 一个序列的DFT，实际上就是这个序列的频谱（DTFT）在一个周期内等间隔抽样的样点值。
 通过这种方法，我们可以通过DFT来观察序列频谱（DTFT）在一个周期内的样点，从而分析信号的频谱。
 DFT和IDFT实际上就是DFS在一个周期上的取值。
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
 
 
5.系统函数，频率响应
 
幅频特性、相频特性、差分方程
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
重点解题知识
 

 
 
 
 
 
 
典型题
 

 
 
 
 
 
离散LTI系统的系统函数
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
6.数字信号系统流图直接型画法（符号+、-）
 

 
 
 
 
 
 
 
 
 
 
 

 
 
 

 
 
 
 
 
 
 
7.稳定、因果、左边、右边、双边
 
因果、稳定：
 通过求单位脉冲响应h[k]
 因果性：判断h[k]=0，k<0是否成立
 稳定性：
 
 
 
 化简难题：
 
 左边、右边、双边：
 
 
 
 
 
 
 
 
8.z变换求解 z的反变换（部分分式，留数法）
 
h(n) H(z) x(t)/z x(z)z^n-1
 
 
 
 
 
 
 
 
 
 
 
9.时间分辨率，频率分辨率，信号最高截止频率
 
To、N、w、f、N、k 量纲
 
 
 
 
 
 
 
 

 
基于MATLAB的语音信号处理
 


 
摘要：语音信号处理是目前发展最为迅速的信息科学研究领域中的一个，是目前极为活跃和热门的研究领域，其研究成果具有重要的学术及应用价值。语音信号处理的研究，对于机器语言、语音识别、语音合成等领域都具有很大的意义。MATLAB软件以其强大的运算能力可以很好的完成对语音信号的处理。通过MATLAB可以对数字化的语音信号进行时频域分析，方便地展现语音信号的时域及频域曲线，并且根据语音的特性对语音进行分析。本文主要研究了基于MATLAB软件对语音信号进行的一系列特性分析及处理，帮助我们更好地发展语音编码、语音识别、语音合成等技术。本文通过应用MATLAB对语音信号进行处理仿真，包括短时能量分析、短时自相关分析等特性分析，以及语音合成等。
 
关键词：语音信号；MATLAB；特性分析；语音合成

引言
        人类交换信息最方便的、最快捷的一种方式是语言。在高度发达的信息社会中，用数字化的方法进行语音的识别、合成、增强、传送和储存等是整个数字化通信网中最重要、最基本的组成部分之一。数字电话通信、高音质的窄带语音通信系统、智能机器人、声控打字机、语言学习机、自动翻译机等，都要用到语音信号处理技术，随着现在集成电路和微电子技术的飞速发展，语音信号处理系统逐步走向实用化[1]。
        语音信号处理是一个新兴的交叉学科，是语音和数字信号处理两个学科的结合产物。与认知科学、心理学、语言学、计算机科学、模式识别和人工智能学科有着密切的联系。语音信号处理技术的发展依赖于这些学科的发展，语音信号处理技术的进步也将促进这些领域的进展。语音信号处理目的是得到一些语音特征参数，以便高效的传输或存储，或通过某种处理以达到特定目的，如语音合成，辨识出讲话者、识别出讲话的内容等。随着现代科学技术和计算机技术的发展，除了人与人的自然语言的沟通，人机对话和智能机领域也开始使用语言。这些人造的语言拥有词汇，语法，语法结构和语义内容等。
        语音信号处理的研究可以追溯到1876年贝尔电话的发明，其在真正意义上首次用声电，电声转换技术实现了远距离语音传输。 1939年Homer Dudley提出并研制成功第一个声码器，奠定了语音产生模型的基础，其在语音信号处理领域具有划时代的意义。在20世纪40年代，一种语言声学的专用仪器语谱图仪问世。它可以让你把语音的时变频谱用语图表示出来，得到一个“可见的语言”。 1984年哈斯金斯实验室研制成功语音回放机，此仪器可以自动转换手工绘制的语谱图成为语言，并进行语音合成。随着计算机的出现，语音分析技术可以在计算机上进行。此时语音信号处理无论是在基础研究或在技术应用，都已取得了突破性进展。现在语音信号可分为三个主要分支，即语音编码，语音识别和语音合成技术[10]。
        语音编码技术。语音编码的目的就是在保证一定语音质量的前提下，尽可能降低编码比特率来节省频率资源。语音编码技术的研究开始于1939年， Homer Dudley提出并实现了在低带宽电话电报上传输语音信号的通道声码器，第一个研制成功了语音编码技术。到20世纪70年代，国际电联于1972年发布了64kbit/s脉冲编码调制（PCM）语音编码算法的G.711建议，它被广泛应用于数字交换机、数字通信等领域，从而占据统治地位。在1995年11月ITU-T SG15全会上共轭代数码激励线性预测（CS-ACELP）的8kbit/s语音编码G.729建议被通过，并于1996年6月ITU-T SG15会议上通过G.729附件A：减少复杂度的8kbit/s CS-ACELP语音编解码器，正式成为国际标准[1]。
        语音识别技术。语音识别的研究开始于20世纪50年代贝尔实验室的Audry系统，它是第一个可以识别10个英文数字的语音识别系统， 1959年Fry和Denes等人采用频谱分析和模式匹配来进行识别决策构建音素识别器来辨别9个辅音和4个元音。20世纪60年代末单语音识别的研究取得实质性进展，并将其作为一个重要的课题。一方面是因为计算机的计算能力有了迅速的提高，计算机能够提供实现复杂算法的硬件、软件；另一方面，数字信号处理在当时有了蓬勃发展，从而自20世纪60年代末开始引起了语音识别的研究热潮。
        语音合成技术。第一个合成器是在1835年由W.von Kempelen发明，经过Weston改进的机械讲话机。机器完全模仿人的生理过程，分别应用了特别设计的哨和软管模拟肺部空气动力和口腔。Homer Dudley在1939年发明了第一台电子语音合成器，它不是一个简单的生理过程的模拟，而是在电子电路基础上来实现语音产生源。本文关于语音信号处理方面主要研究了语音合成。语音合成已经在许多方面得到了实际应用，方便了人们的生活，创造了很好的经济效益和社会效益，如公共交通中的自动报站、各种场合的自动报警、电话自动查询服务、文本校对中的语音提示等。综观语言合成技术的研究，语音合成发展方向为提高合成语音的自然度、丰富合成语音的表现力、降低语音合成技术的复杂度等。
一、语音信号处理基本知识与仿真环境介绍
1.1 语音信号处理基本知识
1.1.1语音信号分析技术
        语音信号分析是语音信号处理的前提和基础，只有分析出可表示语音信号本质特征的参数，才有可能利用这些参数进行高效的语音通信、语音合成和语音识别等处理。而且，语音合成的音质好坏，语音识别率的高低，也都取决于对语音信号分析的准确性和精确性。因此语音信号分析在语音信号处理应用中具有举足轻重的地位。
        贯穿于语音分析全过程的是“短时分析技术”。语音信号从整体来看其本质特征的参数是随时间而变化的，所以它是一个非稳态过程，不能用处理稳态信号的数字信号处理技术对其进行分析处理。但是，由于不同的语音是由人的口腔肌肉运动构成声道某种形状而产生的响应，而这种口腔肌肉运动相对于语音频率来说是非常缓慢的，所以从另一方面看，虽然语音倍号具有时变特性，但是在一个短时间范围内(一般认为在10～30ms的短时间内)，其特性基本保持不变即相对稳定，因此可以将其看作是一个准稳态过程，即语音信号具有短时平稳性。所以任何语音信号的分析和处理必须建立在“短时”的基础上．即进行“短时分析”，将语音信号分为一段一段来分析其特征参数，其中每一段称为一“帧”，帧长一般取为10～30ms。这样，对于整体的语音信号来讲，分析出的是由每一帧特征参数组成的特征参数时间序列[4]。
        根据所分析参数的性质的不同，可将语音信号分析分为时域分析、频域分析、倒领域分析等；时域分析方法具有简单、计算量小、物理意义明确等优点，但由于语音信号最重要的感知特性反映在功率谱中，而相位变化只起着很小的作用，所以相对于时域分析来说频域分析更为重要。
1.1.2语音信号处理理论依据
    采样定理。在进行模拟/数字信号的转换过程中，当采样频率大于信号中最高频率的2倍时，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5~10倍。采样定理又称奈奎斯特定理。
        采样位数。采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。
    采样频率。样频率是指计算机每秒钟采样多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据奈奎斯特理论，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。
        采样位数与采样频率的关系。采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度，每增加一个采样位数相当于力度范围增加了6dB，采样位数越多则捕捉到的信号越精确，对于采样率来说你可以想象它类似于一个照相机，44.1khz意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样率越高，计算机提取的图片越多，对于原始音频的还原也越加精确。
1.2 实现平台MATLAB 7.0介绍
1.2.1 MatLab软件基本介绍
        MATLAB产生于1982年，是一种效率高、功能强的数值计算和可视化计算机高级语言，它将信号处理、数值分析和图形显示结合一体，形成了一个极其方便又强大的操作环境，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言的编辑模式，代表了当今国际科学计算软件的先进水平[7]。
        MATLAB7.0是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。MATLAB 7.0 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。
1.2.2 MatLab与语音处理的关系
        MATLAB软件以其强大的运算能力可以很好的完成对语音信号的处理。通过MATLAB可以对数字化的语音信号进行时频域分析，方便地展现语音信号的时域及频域曲线，并且根据语音的特性对语音进行分析。例如，请浊音的幅度差别、语音信号的端点、信号在频域中的共振峰频率、加不同窗和不同窗长对信号的影响、LPC分析、频谱分析等[3]。
        同时，通过MATLAB可以对数字化的语音信号进行估计和判别。例如，根据语音信号的短时参数，一级不同语音信号的短时参数的性质对一段给定的信号进行有无声和请浊音的判断、对语音信号的基音周期进行估计等。另外，通过利用MATLAB编程可以对语音信号进行各种处理。由于MATLAB是一种面向科学和工程计算的高级语言，允许用数学形式的语言编程，又有大量的库函数，所以编程简单、编程效率高、易学易懂，我们可以对信号进行加噪去噪、滤波、截取语音等，也可进行语音编码、语音识别、语音合成的处理等。总之，对于语音信号进行分析处理，MATLAB软件是当今比较高效的平台。
二、语音信号的特点与采集
2.1语音信号的特点分析
        语音信号的特点可以分为时域方面和频域方面。
        在时域内，语音信号具有“短时性”的特点，即在总体上，语音信号的特征是随着时间而变化的，但在一段较短的时间间隔内，语音信号保持平稳。
        在频域内，语音信号的频谱分量主要集中在300～3400Hz的范围内。利用这个特点，可以按8kHz的采样率对语音信号进行采样，得到离散的语音信号。语音信号的这两种特点均可通过MATLAB软件表现出来，如图2.1和图2.2所示。


图2.1 语音信号时域图


图2.2 语音信号频域分析
2.2语音信号的采集
2.2.1语音信号的量化编码采样
        在将语音信号进行数字化前，必须先进行防混叠预滤波，预滤波的目的有两个，一是抑制输入信导各领域分量中频率超出/2的所有分量(为采样频率)，以防止混叠干扰；二是抑制50Hz的电源工频干扰。这样，预滤波器必须是一个带通滤波器，设其上、下截止颜率分别是和，则对于绝人多数语音编译码器，=3400Hz、＝60~100Hz、采样率为＝8kHz；而对于语音识别而言，当用于电话用户时，指标与语音编译码器相同。当使用要求较高或很高的场合时＝4500Hz或8000Hz、＝60Hz、＝10kHz或20kHz。
        为了将原始模拟语音信号变为数字信号，必须经过采样和量化两个步骤，从而得到时间和幅度上均为离散的数字语音信号。采样也称抽样，是信号在时间上的离散化，即按照一定时间间隔△t在模拟信号x(t)上逐点采取其瞬时值。采样时必须要注意满足奈奎斯特定理，即采样频率必须以高于受测信号的最高频率两倍以上的速度进行取样，才能正确地重建波它是通过采样脉冲和模拟信号相乘来实现的。
        在采样的过程中应注意采样间隔的选择和信号混淆：对模拟信号采样首先要确定采样间隔。如何合理选择△t涉及到许多需要考虑的技术因素。一般而言，采样频率越高，采样点数就越密，所得离散信号就越逼近于原信号。但过高的采样频率并不可取，对固定长度（T）的信号，采集到过大的数据量（N=T/△t），给计算机增加不必要的计算工作量和存储空间；若数据量（N）限定，则采样时间过短，会导致一些数据信息被排斥在外。采样频率过低，采样点间隔过远，则离散信号不足以反映原有信号波形特征，无法使信号复原，造成信号混淆。根据采样定理，当采样频率大于信号的两倍带宽时，采样过程不会丢失信息，利用理想滤波器可从采样信号中不失真地重构原始信号波形。量化是对幅值进行离散化，即将振动幅值用二进制量化电平来表示。量化电平按级数变化，实际的振动值是连续的物理量。具体振值用舍入法归到靠近的量化电平上。
        语音信号经过预滤波和采样后，由A/D变换器变换为二进制数字码。这种防混叠滤波通常与模数转换器做在一个集成块内，因此目前来说，语音信号的数字化的质量还是有保证的。市面上购买到的普通声卡在这方面做的都很好，语音声波通过话筒输入到声卡后直接获得的是经过防混叠滤波、A/D变换、量化处理的离散的数字信号。
2.2.2利用Windows录音器采集语音信号
在本次设计中，可以利用Windows自带的录音机录制语音文件，图2.3是基于PC机的语音信号采集过程，声卡可以完成语音波形的A/D转换，获得WAV文件，为后续的处理储备原材料。调节录音机保存界面的“更改”选项，可以存储各种格式的WAV文件。
  

图2.3 基于PC机的语音采集过程
第三章 语音信号的分析
3.1 语音信号的短时能量分析

        一定时宽的语音信号，其能量的大小随时间有明显的变化。清音信号和浊音信号之间的能量差别相当显著。其中清音段（以清音为主要成份的语音段），其能量比浊音段小得多[10]。因此，对语音的短时能量进行分析，可以描述语音的这种特征变化情况。定义短时能量为如式（3-1）所示。

                                                
（3-1）

其中N为窗长。特殊地，当采用矩形窗时，可简化为如式（3-2）所示。

                                                        
（3-2）

也可以从另外一个角度来解释。令

                                                          
（3-3）                                         

则 可表示为如式（3-4）所示。

                                          
（3-4）  

        可以理解为，首先语音信号各个样点值平方，然后通过一个冲击响应为h(n)的滤波器，输出为由短时能量构成的时间序列。

        短时能量的计算直接受冲击响应的选择即窗函数的选择的影响。如果冲击响应的幅度是恒定的，它的序列长度N（即窗长）会很长，将其等效为非常窄的低通滤波器，这时冲击响应对 产生的平滑的作用比较明显，使短时能量基本没有很大的变化，将不能表现出语音的时变的特性。相反，如果冲击响应的序列长度过于小，等效窗就不能提供出够用的平滑，以导致语音的振幅在瞬时的变化的许多细节仍被留了下来，进而不能看出振幅包络变化的规律，一般我们要求窗长是几个基音周期的数量级。

        图3.1为采样率8000kHZ，16位，单声道的一个语音信号（单词“earth”）在不同矩形窗长时的短时能量函数，我们会发现：语音信号的幅度变化在被短时能量所反映时，窗长的长短都有影响。



  
 


 图3.1 不同矩形窗长的短时能量函数
        我们知道，单词earth前半部分是浊音，后半部分是清音。由以上分析结果可知，浊音部分的能量较之清音部分要大得多，而清音部分的能量相当小，几乎为零。

        对语音信号进行短时能量函数运算，可实现以下三点应用：

（1）可用于区分清音段与浊音段。En值大对应于浊音段，En值小对应于清音段。

（2）可用于区分浊音变为清音或清音变为浊音的时间（根据En值的变化趋势）。

（3）对高信噪比的语音信号，也可以用来区分有无语音（语音信号的开始点或终

止点）。无信号（或仅有噪声能量）时，En值很小，有语音信号时，能量显著增大。


3.2短时自相关分析

        对于确定性信号序列，自相关函数定义如式（3-5）所示。

                                                  
（3-5）

        对于随机性信号序列或周期性信号序列，自相关函数的定义如式（3-6）所示。

                                            
（3-6）

        自相关函数具有以下几项性质：

    (1)若序列是周期性的，假设序列周期为 ，那么其自相关函数也是具有相同周期的周期函数，即 

    (2)自相关函数是偶函数，即R(k)=R(-k)；

    (3)当k=0时，自相关函数有极大值，即 
；

    (4)R(0)为随机性序列的平均功率或确定性信号序列的能量。

        自相关函数的上述性质，完全可以适用于语音信号的时域分析中。例如，浊音语音波形序列具有周期性，因此可用自相关函数求出这个周期，即是基音周期。此外，自相关函数也可用在语音信号的线性预测分析中。

短时自相关函数的定义如式（3-7）所示。 

                                   
（3-7）

令 ，并且 ，可以得到如下式子，如（3-8）所示。

                             
（3-8）

        如图3.2是在不同的矩形窗窗长条件下单词earth的语音自相关的函数的波形。

        对两图分析可得：清音信号的短时自相关函数的波形不具有周期性，也没有明显的峰值，且随着延时k的增大迅速变小，因此其接近于随机噪声；浊音是具有周期性的信号，浊音信号的周期为自相关函数的周期，由此可知，语音信号的性质是浊音还是清音，如果是浊音，还可以得出它的基音周期，它的基音周期可由自相关函数波形中的第一个峰值的位置来估计。所以，自相关函数常用作一下两种作用：

(1) 区分语音信号是清音还是浊音；

(2) 估计浊音语音信号的基音周期[4]。



  
 

图3.2 不同的矩形窗窗长下短时自相关

第四章 语音合成
4.1 语音合成技术概述

4.1.1 语音合成技术的意义                                                   

        语音合成技术涉及声学、语言学、数字信号处理技术、多媒体技术等多个领域, 是当今世界强国竞相研究的热门技术之一。语音合成技术可分为参数合成和波形拼接两种方法。早期的研究主要是采用参数合成方法, 它是计算发音器官的参数, 从而对人的发音进行直接模拟。语音合成已经在许多方面得到了实际应用，方便了人们的生活，创造了很好的经济效益和社会效益，如公共交通中的自动报站、各种场合的自动报警、电话自动查询服务、文本校对中的语音提示等[8]。

        本文主要利用载波调制技术进行语音合成。基于载波调制的语音信号合成是以语音信号处理技术、数字信号处理技术为基础，依托于电子计算机、Windows操作系统、MATLAB处理软件等工具将两个信号合成为一个信号。具有较强的实用性、可操作性等特点。

4.1.2 基于载波调制语音合成的基本原理

            语音信号合成是一个“分析—存储—合成”的过程。一般是选择合适的基本单元，将基本单元用一定的参数编码方式或波形方式进行存储，形成一个语音库。合成时，根据待合成语音信号，从语音库中取出基本单元进行合成，并将其还原成语音信号。在语音合成中，为了便于存储和后续分析，必须先将语音信号进行预分析、预处理、波形变换等一系列操作。其中，基元是语音合成处理的最小单元，待合成的语音库就是所有语音基元的某中集合。根据基元的选择方式以及其存储形式的不同，可以将合成方法笼统的分为波形合成方法和参数合成方法。

        波形合成是一种相对简单的语音合成技术。它把人的发音波形直接存储或者进行进行简单的波形编码后存储，组成一个合成的语音库；合成时，根据待合成的信息，在语音库中取出相应单元的波形数据，拼接或编辑到一起，经过解码还原成语音。该语音合成技术具有一定的局限和不足，但对语音信号具有数据量庞大的特点，这种误差在某种范围内是可以接受的。

        基于载波调制的语音信号合成是基于信号的振幅调制原理而实现的。将低频信号加载到高频载波信号的过程，或者说把信息加载到信息载体上以便传输的处理过程，称为调制。所谓“加载”，其实质是使高频载波信号（信息载体）的某个特性参数随信息信号幅值的大小程线性变化。基于载波调制的语音信号合成是以语音一信号作为调制信号，语音二信号为载波信号来进行合成一种以语音二信号声色表述语音一内容的新信号。这种调制方式是用传递的低频信号（如代表语言、音乐、图像的电信号）去控制作为传送载体的高频振荡波（称为载波）的幅度，是已调波的幅度随调制信号的大小线性变化，而保持载波的角频率不变。

4.2 基于载波调制的语音合成基本知识

4.2.1 关键函数希尔伯特变换介绍

        本文语音合成的设计思路是用一个语音信号的包络去调制另一个语音信号，实现语音的合成。这就用到了一个关键的函数，希尔伯特变换。在数学与信号处理的领域中，一个实值函数的希尔伯特变换是将信号s(t)与1/(πt)做卷积，以得到s'(t)。因此，希尔伯特变换结果s'(t)可以被解读为输入是s(t)的线性时不变系统的输出，而此系统的脉冲响应为1/(πt)。这是一项有用的数学工具，用在描述一个以实数值载波做调制的信号之复数包络，出现在通讯理论中发挥着重要作用[9]。

        希尔伯特变换的频域数学表达式如式（4-1）所示。                       

                                                      
（4-1）

其中F是傅里叶变换，i是虚数单位，ω是角频率。

        希尔伯特变换等效于 π/2的相移，对正频率产生-π/2的相移，对负频率产生π/2相移，或者说，在时域信号每一频率成分移位1/4波长，因此，希尔伯特变换又称为90度移相器。

        MATLAB提供了计算Hilbert变换的函数，其格式为y=Hilbert(x)。但需注意的是，该函数计算出的结果是序列的解析信号，其虚部才是序列的Hilbert变换。

        希尔伯特变换在语音信号处理中具有两个性质： 序列x(n)通过Hilbert变换器后，信号频谱的幅度不发生变化，这是因为Hilbert变换器是全通滤波器，引起频谱变化的只是其相位； 序列x(n)与其Hilbert变换是正交的[6]。

4.2.2 信号调制

        所谓调制，就是将调制信号加载在三个参数中的某一个参数上，或幅值、或频率、或相位，随调制信号大小成线性或非线性变化的过程。主要有三种基本调制方法，第一种是把调制信号加载在载波信号的幅值上，称为幅度调制 ，简称AM；第二种是把调制信号加载在载波的频率上，称为频率调制，简称FM。 第三种是把调制信号装载在载波的相位上，称为相位调制，简称PM[10]。 本设计采用的是第一种方法，用采集到的语音二信号去对语音一信号进行幅度调制，实现语音合成的目的。

        采用调幅调制是因为其以下特点在语音信号处理中得到很好的应用。一是调幅波的振幅（包络）随调制信号变化，而且包络的变化规律与调制信号波形一致，表明调制信号（信息）记载在调幅波的包络中；二是调制系数反应了调幅的强弱程度，一般情况下，调制系数越大调幅度越深。

        当调制系数为0时，表示未调幅，即无调幅作用；

        当调制系数为1时，此时包络的振幅最小值为0；

        当调制系数大于1时，已调波的包络与调制信号不一样，产生严重的包络失真，称为过量调幅。

4.3 语音信号合成过程

4.3.1 语音信号合成流程图

        用MATLAB 处理音频信号的基本流程是：先将WAV 格式音频信号经wavread 函数转换MATLAB 列数组变量；再用MATLAB 强大的运算能力进行数据分析和处理，如时域分析、频域分析、数字滤波、信号合成、信号变换、识别和增强等等；处理后的数据如是音频数据，则可用wavwrite函数 转换成WAV 格式文件或用sound、wavplay 等函数直接回放。本设计的语音合成流程图如图4.1所示。



  
 


 图4.1  语音信号合成流程图
4.3.2 语音信号的采集

        分析和处理音频信号，首先要对声音信号进行采集，MATLAB 的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可直接控制声卡进行数据采集。Windows 自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WAV 格式文件，供MATLAB 相关函数直接读取、写入或播放。本文以WAV 格式音频信号作为分析处理的输入数据。

4.3.3 语音信号的合成

        声音信号是一种非平稳信号，如果采用简单的时变系统的分析方法，将会产生很大的失真，但是在一小段时间内声音信号完全可以视为平稳信号。因此必须对语音信号做预处理。在本次语音信号合成中采用加窗截断，分帧处理将非平稳信号近似转换为平稳信号。

        声音信号特征量提取。声音信号特征量提取包括语音一声音信号声色（频率）的提取和语音二声音信号包络的提取。语音二声音信号包络的提取采用希尔伯特变换实现，得到语音二声音信号的复数包络。

        获取语音信号起始位置。在录音过程中控制两段声音从相同的起点开始录取并不是一件容易的事，但是如果不确定语音信号的起始位置直接对语音信号进行合成既存在数据量大又会带来较大的误差。本设计语音合成中拟定连续四个时间点的语音信号强度不为0，则认为语音信号开始，也即找到信号起始位置。

        语音信号合成。语音信号合成即是一个调幅载波的过程，是以语音一信号作为载波信号，语音二信号包络作为调幅信号来实现语音合成。实际的载波是一个物理可实现的复杂过程，本语音合成中采用语音一信号点乘包络信号，实际的载波是一个物理可实现的复杂过程，并非简单地乘积运算，然而，得到的合成声音信号效果并不理想，但其波形仍能反应载波过程的实质。

4.4 语音信号合成结果及分析

4.4.1 语音信号预处理结果及分析

        该处理过程以语音一信号和语音二声音信号为分析样本。使用Windows系统自带录音器分别录下语音一和语音二，分别命名为one和two,保存为WAV格式。通过MATLAB对所录语音进行采样，采样频率 为16000Hz，获取语音信号并进行加窗。语音一和语音二的时域波形图如图4.2所示，时域图反映出了语音信号的非平稳性。

        对采集到的语音信号分别做傅里叶变换进行频谱分析，并显示频谱图，观察各自的幅频谱特性。语音一和语音二的声音信号幅频特性如图4.3所示，语音一和语音二的声音FFT图如图4.2和图4.3。该频谱图横坐标并未进行对应关系处理，但仍不失其频谱特性的本质，由频谱图可清楚地看到样本声音主要以低频为主。人的语音信号频率一般集中在1kHz之前，从声音频谱的包络来看, 根据采样定理，信号宽度近似取为1kHz，重放语音后仍可较清晰的听出原声, 不存在声音混叠现象。



  
 


 图4.2  信号预处理之后时域图


  
 


 图4.3 信号预处理之后频域图
4.4.2 合成语音信号结果及分析

        合成语音信号的实质是用语音二信号的包络调制语音一信号振幅的结果。语音二信号包络提取结果如图4.4，该图是语音二信号经希尔伯特—黄变换的虚部显示，因为希尔伯特—黄变换是一个时域信号与1/（πt）的卷积，其结果是载波做调制信号之复数包络，必然蕴含虚部成分，取其虚部的结果必然与时域信号有着直观上的差别，但仍是信号的包络成分。



  
 


 图4.4 语音二信号包络图
        合成信号的时域显示结果如图4.5所示，该合成信号是以语音一信号的特性和语音二信号的幅度变化的，由其快速傅里叶变换的结果更证实了这一点，其幅频特性与语音一信号的幅频特性更接近。



  
 


 图4.5 合成语音信号的时域波形


  
 


 图4.6 合成语音信号的幅频特性



  
 

图4.7 合成语音信号快速傅里叶变换结果


结  论
        随着语音技术的逐渐成熟,语音信号处理技术也在不断发展,不断完善。本文主要研究了通过对语音信号短时能量、短时自相关等特性参数的分析，使我进一步了解了语音信号的特性，明白了只有准确分析并提取出语音信号的特征参数，才能够利用这些参数进行语音编码、语音合成等处理。另外在语音处理方面，我选择了语音合成这一处理方式。基于载波调制的语音处理实现简单，运用广泛，研究这一语音合成方法及特性，对于更加深入地进行各种语音处理有着重要的意义。这次设计我是通过了MATLAB这一平台，MATLAB软件以其强大的运算能力可以很好的完成对语音信号的处理，因此，近一步的加强对MATLAB的研究对我以后的学习会起到很大的帮助。

        至此，设计基本符合要求。但是由于个人能力的有限，采集的语音信号清、浊音区分不明显，导致对语音进行短时自相关得出的波形特征不明显。考虑解决方案是通过专业的设备采集语音信号。除此之外，本设计必有其他欠妥之处，请各位老师给予指正！




参考文献
[1] 张雪英.数字语音处理[M].北京:电子工业出版社, 2010.

[2] 郑君里,应启绗,杨为理.信号与系统[M].北京:高等教育出版社,2000.

[3] 薛年喜.MATLAB在数字信号处理中的应用[M].北京:清华大学出版社,2003.

[4] 胡航.语音信号处理[M].哈尔滨:哈尔滨工业大学出版社, 2000.

[5] 易克初,田斌,付强.语音信号处理[M].北京:国防工业出版社, 2000.

[6] 万永革.数字信号处理的MATLAB实现[M].北京:科学出版社,2007.

[7] 刘卫国. MATLAB程序设计与应用[M].北京:高等教育出版社,2006.

[8]王嘉梅.基于MATLAB的数字信号处理与时间开发[J].西安:西安电子科技大学出版社,2007:10-14.

[9] 程佩青.数字信号处理教程（第二版）[M].北京:清华大学出版社,2010.

[10] 韩纪庆 张磊 郑铁然.语音信号处理[M].北京:清华大学出版社,2004.

[11] 徐明远,邵玉斌. Matlab仿真在通信与电子工程中的应用[M].西安:西安电子科技大学出版社,2005. 

[12] 邓华. Matlab通信仿真及应用实例详解[M].北京:人民邮电出版社,2005. 

[13] 张照明,刘政波,刘斌等. 应用Matlab实现信号分析处理[C].北京:科学出版社,2006.

[14] 徐守时. 信号与系统理论方法和应用[M].合肥:中国科学技术大学出版,1999. 


[15] 高俊斌. Matlab语言与程序设计[M].武汉:华中理工大学出版社,1998. 
附  录

附录A   语音信号特性分析程序
%语音信号时域频域显示%
[y,Fs,bits]=wavread('biye.wav');%读出信号、采样率和采样位数
y=y(:,1);%取单声道
sigLength=length(y);
Y=fft(y,sigLength); 
Pyy=Y.* conj(Y) / sigLength;
halflength=floor(sigLength/2);
f=Fs*(0:halflength)/sigLength;
figure;plot(f,Pyy(1:halflength+1));
xlabel('Frequency(Hz)');
t=(0:sigLength-1)/Fs;
figure;
plot(t,y);
xlabel('Time(s)');

%语音信号短时能量%
x=wavread('biye.wav');
%x=fscanf(fid,'% f');
%fclose(fid);
s=fra(50,25,x)
s2=s.^2;
energy=sum(s2,2)
subplot(2,2,1)
plot(energy);
xlabel('帧数')
ylabel('短时能量 E')
legend('N=50')
%axis({0,1500,0,10*10^5})
s=fra(100,50,x)
s2=s.^2;
energy=sum(s2,2)
subplot(2,2,2)
plot(energy);
xlabel('帧数')
ylabel('短时能量 E')
legend('N=100')
%axis({0,750,0,2*10^6}) 
s=fra(400,200,x)
s2=s.^2;
energy=sum(s2,2)
subplot(2,2,3)
plot(energy);
xlabel('帧数')
ylabel('短时能量 E')
legend('N=400')
%axis({0,190,0,7*10^6})
s=fra(800,400,x)
s2=s.^2;
energy=sum(s2,2)
subplot(2,2,4)
plot(energy);
xlabel('帧数')
ylabel('短时能量 E')
legend('N=800') 
%axis({0,95,0,14*10^6})

%语音信号短时自相关%
x=wavread('biye.wav');
s1=x(1:320);
N=320;   %选择的窗长，加N=320的矩形窗
A=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m)*s1(m+k-1);   %计算自相关
end
A(k)=sum;
  end
for k=1:320
A1(k)=A(k)/A(1);        %归一化A（k）
   end
N=160;                  %选择的窗长，%加N=160的矩形窗
B=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m+k-1);   %计算自相关
end
B(k)=sum;
end
for k=1:320
B1(k)=B(k)/B(1);      %归一化B（k）
end
N=70;                 %选择的窗长，加N=70的矩形窗
C=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m)*s1(m+k-1);        %计算自相关
end
C(k)=sum;
end
for k=1:320
C1(k)=C(k)/C(1);                %归一化C（k）
end
s2=s1/max(s1)
figure(1)
subplot(4,1,1)
plot(s2)
title('语音信号')
xlabel('样点数')
ylabel('幅值')
axis([0,320,-2,2])
subplot(4,1,2)
plot(A1)
xlabel('延时k')
ylabel('R(k)')
axis([1,320,-2,2]);
legend('N=320')
subplot(4,1,3)
plot(B1);
xlabel('延时k')
ylabel('R(k)')
axis([1,320,-2,2]);
legend('N=160')
subplot(4,1,4)
plot(C1);
xlabel('延时k')
ylabel('R(k)')
axis([0,320,-2,2]);
legend('N=70')
附录B  语音合成主程序
[y1,fs,bits]=wavread('one');      %读取语音一信号
[y2,fs,bits]=wavread('two');      %读取语音二信号
L1=length(y1);                    %测定语音一信号长度
L2=length(y2);                    %测定语音二信号长度
a1=y1.*hamming(L1);               %加窗预处理
a2=y2.*hamming(L2);               %加窗预处理
L1=length(a1);                    %测定语音一信号长度
L2=length(a2);                    %测定语音二信号长度
%采样信号的时域显示
figure(1);
subplot(211);
plot(a1);
title('语音一载波信号时域波形');
subplot(212);
plot(a2);
title('语音二调幅信号时域波形');
%傅里叶频谱绘制
F1=fft(a1,L1);                    
F2=fft(a2,L2);
AF1=abs(F1);
AF2=abs(F2);
figure(2);
subplot(211);
plot(AF1);
title('语音一载波信号幅频特性显示');
subplot(212);
plot(AF2);
title('语音二调幅信号幅频特性显示');
figure(3);
freqz(F1);
title('语音一载波信号FFT频谱显示');
figure(4);
freqz(F2);
title('语音二载波信号FFT频谱显示');
%获取语音一信号的开始位置
for i=1:L1-4
     g(i)=a1(i).*a1(i+1).*a1(i+2).*a1(i+3).*a1(i+4);%认为连续4个幅值不为0的信号即为开始
     if g(i)~=0
         break;
     else i=i+1;
     end
end
I=i;
 
% 获取语音二信号开始位置
for j=1:L2-4
     m(j)=a2(j).*a2(j+1).*a2(j+2).*a2(j+3).*a2(j+4);
     if m(j)~=0
         break;
     else j=j+1;
     end
end
J=j;
%语音二信号hilbert变换
H=hilbert(a2);
figure(5);
plot(abs(H));
title('语音二信号包络显示');
%信号对齐,语音二包络调制语音一振幅
max1=max(I,J);
for k=1:L1-max1
    N(k)=a1(i).*H(j);
    i=i+1;
    j=j+1;
end
%N=N';
N = N/(max(abs(N)) * 1.05);
wavwrite(N,16000,16,'HC.wav');
figure(6);
plot(imag(N));
title('合成信号时域显示');
pause(1);
sound(10*N,fs);
FN=fft(N);
figure(7);
freqz(FN);
title('合成声音信号FFT显示');
figure(8);
plot(abs(FN));
title('合成声音信号的幅频特性');

~~ 如果有什么问题可以在我的个人博客留言 ，我会及时回复。欢迎来访交流 ~~
 

               IIR数字滤波器设计（数字信号处理）
 
一、实验目的
 
1.熟悉双线性变换法设计IIR数字滤波器的原理与方法。
 
2.掌握IIR数字滤波器的MATLAB实现方法，会调用ellipord()和ellip() 
 
函数设计各种滤波器。
 
3.观察分析滤波器输入输出数据波形，理解数字滤波的概念。
 
             
 
二、实验原理及步骤
 
（一）实验原理-双线性变换法
 
数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性:
 
 
其中、 分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性），是数字滤波器的单位取样响应的频谱，又称为数字滤经过滤波后的频域响应。只要按照输入信号频谱的特点和处理信号的目的，适当选择，使得滤波后满足设计的要求，这就是数字滤波器的滤波原理。
 
数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为：
 
 
系统函数为：
 
             
 
设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。
 
基本设计过程如下：
 
1.将给定的数字滤波器指标转换成模拟滤波器的指标。
 
2.设计模拟滤波器。
 
3.将模拟滤波器转换成数字滤波器系统函数。
 
（二）实验步骤
 
1.在完成滤波器设计之后，采用filter()对输入信号进行滤波处理。调用信号产生函数mstg产生由抑制载波调制信号相加构成的复合信号st，该函数还会自动绘图显示其时域波形和幅频特性曲线，如下图1所示。
 
 
 
由图1中（a）和（b）可见，三路信号时域混叠无法在时域分离，但在频域是可以分离的，可以通过滤波的方法进行分离，即通过设计IIR滤波器，分离这三个不同频率的信号。
 
2.要求将三路信号进行分离，通过观察st信号的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB。
 
3.调用ellipord()和ellip()分别设计这三个椭圆滤波器，并绘图显示其损耗函数曲线分别如图2，图3，图4。
 
 
                                                                 图2 低通损耗函数曲线
 
 
 
4.调用filter()，用三个滤波器分别对信号产生函数mstg产生的信号进行滤波，分离出st中的不同载波频率信号，并绘图显示三个信号的时域波形，分别如图5，图6，图7。
 
 
 
 
 
 
三、实验结果分析
 
（一）函数mstg
 
阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。第1路调幅信号的载波频率fc1=1000Hz，调制信号频率fm1=100Hz。第2路调幅信号的载波频率fc2=500Hz，调制信号频率fm2=50Hz。第3路调幅信号的载波频率fc2=250Hz，调制信号频率fm2=25Hz。
 
采样点数对频谱图的影响
1.调幅信号
 
当N=800,1600,1800,2000时，调幅信号产生的频谱图如图8,9,10,11所示。
 
 
                                                             图8 N=800时s（t）的波形及频谱图
 
 
 
分析：因为信号s（t）是周期序列，频谱分析时要求观察时间为整数倍周期。s（t）的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25*400Hz，即在25Hz的正弦波的一个周期中采样400个点。所以，当N为400的整数倍时一定为s（t）的整数倍周期。因此，采样点数N=800,1600,2000时，对s（t）进行N点FFT可以得6根理想谱线，而当N=1800时，不是400的整数倍，则不能得到。
 
2.AM调幅信号
 
当N=800,1600,1800,2000时产生的频谱图如图12,13,14,15所示。
 
 
 
 
 
分析：因为信号s（t）时周期序列，频谱分析时要求观察时间为整数倍周期。因此，采样点数N=800,1600,2000时，对s（t）进行N点FFT可以得理想谱线，而当N=1800时，不是400的整数倍，则不能得到。当将该调幅信号修改为AM信号后，s（t）的频谱中有较大的频谱分量。如图所示。
 
IIR滤波器
滤波器参数选取
由（一）可知，三路调幅信号的载波频率分别为250Hz，500Hz，1000Hz。带宽为50Hz，100Hz，200Hz。所以分离混合信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的指标参数选取如下：
 
对载波频率为250Hz的调幅信号，用低通滤波器分离，其通带截止频率fp=280Hz，通带最大衰减ap=0.1dB，阻带截止频率为fs=450Hz，阻带最小衰减as=60dB。对载波频率为500Hz的调幅信号，用低通滤波器分离，其通带截止频率fpl=440Hz、fph=560Hz，通带最大衰减ap=0.1dB，阻带截止频率为fsl=275Hz、fsh=900Hz阻带最小衰减as=60dB。对载波频率为1000Hz的调幅信号，用低通滤波器分离，其通带截止频率fp=890Hz，通带最大衰减ap=0.1dB，阻带截止频率为fs=600Hz，阻带最小衰减as=60dB。
 
 
 
IIR.m
 
function main
st=mstg;   %调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号s
Fs=10000; T=1/Fs;  %采样频率
fp=280; fs=450;wp=2*fp/Fs;ws=2*fs/Fs; %低通
   name='y_1(t)';LHL_filter(wp,ws,name,st,T,'low');
fp_L=440; fp_R=560;fs_L=275; fs_R=900;%带通
   wp=[2*fp_L/Fs,2*fp_R/Fs];ws=[2*fs_L/Fs,2*fs_R/Fs];
   name='y_2(t)'; LHL_filter(wp,ws,name,st,T,'bandpass');
fp=890; fs=600;wp=2*fp/Fs;ws=2*fs/Fs; %高通
   name='y_3(t)';LHL_filter(wp,ws,name,st,T,'high');
end
function LHL_filter(wp,ws,name,st,T,flag)
rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp?
[B,A]=ellip(N,rp,rs,wp,flag); %调用ellip计算椭圆带通DF系统函数系数向量B和A
yt=filter(B,A,st); %滤波器软件实现
 figure; freqz(B,A);%B,A分别为系统函数分子，分母多项式系数向量?
 figure;tplot(st,yt,T,name);
end
function tplot(st,xn,T,name)
%时域序列连续曲线绘图函数?
%?xn:信号数据序列，%?T为采样间隔
N=length(xn);n=0:N-1; t=n*T; Tp=N*T; f=n/Tp;
subplot(311);plot(t,xn);xlabel('t/s');ylabel(name);
axis([0,t(end),min(xn),1.2*max(xn)]);
subplot(312);stem(f,abs(fftshift(fft(st,N))));
title('原 st(t)的频谱');xlabel('f/Hz');ylabel('幅度');
subplot(313);stem(f,abs(fftshift(fft(xn,N))));
title('滤波后xn(t)的频谱');xlabel('f/Hz');ylabel('幅度');
end
function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg 返回三路调幅信号相加形成的混合信号，长度N=1600
N=1600;Fs=10000;T=1/Fs;Tp=N*T; % N为信号st的长度。采样频率Fs=10kHz，Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz                  
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号
st=xt1+xt2+xt3;         %三路调幅信号相加
fxt=fft(st,N);          %计算信号st的频谱
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/2,min(st),max(st)]);title('(a) s(t)的波形');
subplot(212);stem(f,abs(fxt)./max(abs(fxt)));
title('(频谱');xlabel('f/Hz');ylabel('幅度');
end
 
mstg.m
 
function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg; 返回三路调幅信号相加形成的混合信号，长度N=1600
N=2000;Fs=10000;T=1/Fs;Tp=N*T; % N为信号st的长度。采样频率Fs=10kHz，Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz                  
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号
st=xt1+xt2+xt3;         %三路调幅信号相加
fxt=fft(st,N);          %计算信号st的频谱
subplot(2,1,1);plot(t,st);xlabel('t/s');ylabel('s(t)');
title(['N=',num2str(N),'s(t)的波形']);axis([0,Tp/8,min(st),max(st)]);
subplot(2,1,2);stem(f,abs(fxt)/max(abs(fxt)));
title(['N=',num2str(N),'s(t)的频谱']);xlabel('f/Hz');ylabel('幅度');
 axis([0,Fs/5,0,1.2]);
end
 
AM_mstg.m
 
function AM_mstg
N=2000; Fs=10000;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz
xt1=[1+cos(2*pi*fm1*t)].*cos(2*pi*fc1*t);
xt2=[1+cos(2*pi*fm2*t)].*cos(2*pi*fc2*t);
xt3=[1+cos(2*pi*fm3*t)].*cos(2*pi*fc3*t);
st=xt1+xt2+xt3;fxt=fft(st,N);
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title(['N=',num2str(N),' AM s(t)波形图'])
subplot(212);stem(f,abs(fxt)/max(abs(fxt)));title(['N= ',num2str(N),'频谱图'])
axis([0,Fs/5,0,1.2]);xlabel('f/Hz');ylabel('幅度');
end
 
资源点击下载