2022年南邮通达语音时频参数的提取和应用和语音编码实验.pdf

资料收集于网络如有侵权请联系网站删除谢谢精品文档南京邮电大学通达学院实 验 报 告实验名称语音时频参数的提取和应用实验语音编码实验课程名称信息处理技术专业综合实验班级学号姓名开课时间2016/2017 学年, 第 二 学期精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 1 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档实验一语音时频参数的提取和应用一、实验目的熟悉基本的语音时频参数提取方法，观察并比较各类参数采用不同提取方法的效果，了解这些时频参数在语音处理中的应用。二、实验内容1、编写并调试语音时域参数提取程序，包括短时功率、短时过零率等函数。2、编写并调试语音频域参数提取程序，包括频谱、语谱图、基音频率、共振峰等。3、利用提取出的时频参数对语音信号进行静/清 /浊音的划分。三、实验原理语音信号具有短时平稳性，需对其进行分段（帧）处理（帧长一般取10至 30ms之间）。根据激励方式的不同，语音信号可以分成浊音和清音两大类：浊音：激励源可模拟成准周期性的脉冲串。该周期称为基音周期，其倒数就是基音频率，基音频率一般在几百赫兹以内。清音：激励源可模拟成随机白噪声。简化的语音生成数学模型如图所示：准周期脉冲序列发生器随机噪声发生器基音周期时变线性系统浊音清音增益控制时变参数语音取样值对语音时频波形进行绘制与观察，有助于了解语音的清/浊音等特性。对基音周期等声门参数和共振峰等声道参数进行提取与分析，是语音处理与应用的重要内容。对于时变语音信号)(nx，可以计算它的短时功率)(mPx和短时过零率)(mZx。第m帧（帧长为N个样点）的功率计算公式为：21)(1)(mNmnxnxNmP精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 2 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档过零率计算公式为：mNmnxnxsignnxsignNmZ12)1()(1)(对于浊音语音，可以利用其频谱)(X具有丰富的谐波分量的特点，求出其谐波乘积谱，计算公式为：RrrXHPSx1)()(式中， R一般取为5。在谐波乘积谱中，基频分量变得很大，更易于估计基音周期。四、实验方法及程序1. 调用 matlab 中的 wavread 和 wavplay 读入并播放语音文件。2. 调用 fft 命令计算语音的频谱。3. 调用 specgram 命令得到语音的语谱图。4. 调用 plot 命令绘制并观察语音的时频域波形。5. 根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时 功 率 、 短 时 过 零 率 、 基 音 频 率 的 函 数 ： function Px=stpower、 function Zx = stzerocross(x, N) 和function HPSx = hpspectrum(x,N,R)。6. 利用短时功率和短时过零率这两个参数可以对语音信号进行浊音/清音分类。编写实现此分类功能的函数voiunvoi.m。五、实验结果与分析1. 清音和浊音的短时功率、短时过零率各有何特点？清音段的能量一般比浊音段的小得多。浊音时能量集中于较低频率段内，具有较低的过零率，而清音时能量集中于较高频率段内，具有较高的过零率。2. 如何利用短时功率和短时过零率这两个参数对语音信号进行浊音/清音的分类？发浊音时， 由于声门波引起了谱的高频跌落。所以其语音能量集中在3kHz 以下；发清音时，多数能量出现在高频。高频率意味着高的平均过零率，低频率意味着低的平均过零率，一般可认为浊音具有较低的平均过零率，清音具有较高的平均过零率。3. 如何根据时域波形估计各帧元音的基音周期？看图可知四帧语音的周期均为0.01s ，第一帧0100Hz 有 10个峰值点，基频约为10Hz ，第二帧0100Hz 有5个峰值点，基频约为20Hz ，第三帧0100Hz 有 5个峰值点，基频约为20Hz ，第四帧0100Hz 有 4个峰值点，基频约为25Hz 。4. 如何根据基于DFT 的对数幅度谱估计某幀浊音的基音周期？01172 为一个基因周期。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 3 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档5. 如何根据基于DFT 的对数幅度谱，估计出共振峰频率？有 7个峰点，因此基频约为167。6. 时域对语音信号进行加窗，反映在频域，其窗谱对基于DFT 的对数幅度谱有何影响？如何估计出窗谱的主瓣宽度？因 为时 间窗 幅度 的傅 立叶 变换 为脉 冲状 态 ,则频谱偏差 很小 ,当 定常 过程 x(n) 的平均值 为零且 N 很大时ETN()=f()(12)所以加窗的周期图TN( )是频谱密度函数的无偏估计。固有TN()=1+ ( )f()。如果原始信号的频谱成份与FFT 中的谱线完全一致，这种情况下采样数据的长度为信号周期的整数倍，频谱中只有主瓣。7. 如何根据窄带、宽带语谱图，提取出语音的基音和共振峰轨迹？在基音提取中，广泛采用语音波形或误差信号波形的低通滤波。提取共振峰特性最简便的手段是使用语谱仪。提取共振峰还有倒谱法、LPC 分析法等更为有效、准确的方法。（ 做完试验后，记住要回答这些问题）实验过程记录与结果分析1、进入 matlab 1.1 在 Windows 环境中，创建一个名为speech的文件夹，存放语音数据和与实验相关的 Matlab 文件。如：E:speech 。用于实验的语音数据（.mat ）及相应说明（.txt ）包括：digits.mat % digits.txt gliss.mat % 两个包含/i/ 的滑音gliss.txt letters.mat % 英文字母表中26个字母的发音letters.txt ma1.txt ma1-1.mat % 一句连续语音ma1-2.mat % 一句连续语音ma1-3.mat % 一句连续语音ma1-4.mat % 一句连续语音timit.txt timit1.mat % 一句连续语音timit2.mat % 一句连续语音timit3.mat % 一句连续语音timit4.mat % 一句连续语音vowels.mat % 元音 /a/, /i/, /o/, /u/的发音vowels.txt words.mat % 十个孤立字的发音精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 4 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档words.txt 1.2 在 Windows桌面上，双击Matlab的图标，进入Matlab的工作环境。在命令窗口 提示符后，键入cd E:speech指令，进入实验文件夹路径。1.3 在命令窗口键入clear all指令，清除工作空间的所有变量。2、加载语音数据，熟悉matlab 命令2.1 加载“ timit1.mat”语音数据，命令为：load timit1 2.2用命令“who ”和“whos ”列出当前工作空间中的变量。变量“timit1 ”会以包含55911 个元素的矩阵（列向量）形式出现Your variables are: timit1 Name Size Bytes Class Attributes timit1 55911x1 447288 double2.3 利用命令“length ”和“ size”查看语音信号timit1的长度和维数：m = length(timit1) m,n = size(timit1) （粘贴实验结果）m = 55911 m = 55911 n = 1 2.4 通过以下命令，分别查看语音信号timit1中间的1 个、 2 个和 20 个数据值：timit1(5001) timit1(5001,5003) timit1(5001:5020) （粘贴实验结果）ans = 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 5 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档-0.0646 ans = -0.0646 0.0111 ans = -0.0646 0.0307 0.0111 -0.0200 0.2021 0.4012 0.1869 -0.1756 -0.1902 0.0823 0.2363 0.1712 0.0637 -0.0209 -0.0683 -0.0394 0.0081 0.0105 0.0311 0.0696 2.5 Matlab中，符号“ ”可以实现转置运算，运行以下命令进行观察：size (timit1 ) （粘贴实验结果）ans = 1 55911 2.6 使用函数max 和 min 可以分别得到信号的最大值和最小值：max(timit1) min(timit1) （粘贴实验结果）ans = 1 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 6 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档ans = -0.8068 3、语音信号的时域分析3.1 使用以下命令，画出语音信号timit1 （采样频率为16kHz ，长度为3.5 秒）的时域波形。plot (timit1) title ( timit1 ) xlabel ( Time Index,it n ) ylabel ( Amplitude ) axis (0 length(timit1) -1 1) grid 3.2 分别画出一帧浊音和一帧清音的语音时域波形（采样频率为16kHz ，帧长为25ms，每帧有400 个样点） 。subplot(211), plot(timit1(14501:14900) % /a/ in dark.subplot(212), plot(timit1(35501:35900) % /s/ in wash.精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 7 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档3.3 编写个一用以计算信号短时功率的函数“stpower.m ” ，其Matlab代码如下：function Px=stpower (x, N) M = length (x); Px = zeros (M, 1); Px (N)=x (1:N) *x (1:N)/N; For (m = (N+1): M) Px (m) = Px (m-1)+ (x (m) 2-x (m-N) 2)/N; 3.4 编写一个用以计算信号短时过零率的函数“stzerocross.m ” ，其Matlab代码如下：function Zx = stzerocross(x, N) M = length(x); Zx = zeros(M,1); Zx (N+1)=sum (abs (sign (x (2:N+1) sign (x (1:N) ) ) )/ (2*N); for (m= (N+2):M) Zx (m) = Zx (m-1) + (abs (sign (x (m) ) sign (x (m-1) ) ) abs(sign (x (m-N) ) sign (x (m-N-1 ) ) ) )/ (2*N); end 3.5 加载“ digits.mat ”语音数据，该数据包含英文数字“0” 到“ 9”的发音。计算其中的单词“four ”的短时功率和过零率（采样频率为10kHz ，帧长为30ms，每帧有300 个样点） 。load digits; N = 300; x=digits.four1; 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 8 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档Px = stpower (x,N); Zx = stzerocross (x,N); plot(Px*1e -5 Zx x/2000) 3.6 加载“ vowels.mat ”语音数据，该数据包含元音/a/, /i/, /o/, /u/ 的发音。分别画出一帧/a/ 、一帧 /i/ 、一帧 /o/ 和一帧 /u/ 的时域波形（采样频率为10kHz ，帧长为30ms，每帧有300 个样点） 。load vowels subplot(221) plot (vowels.a_1(2001:2300) subplot(222) plot (vowels.i_1(2001:2300) subplot(223) plot (vowels.o_1(2001:2300) subplot(224) plot (vowels.u_1(2001:2300) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 9 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档3.7编写一个利用短时功率和短时过零率这两个参数对语音信号进行浊音/清音分类的函数“voiunvoi.m” ，其Matlab代码如下：function voi = voiunvoi(x,N,Pth,Zth) % Short-time power Px and zero crossing Zx measures. Px = stpower(x,N); Zx = stzerocross(x,N); % Compare estimates with threshold values. voi = (PxPth*max(Px) & (ZxZth); % Shift the voi-flag N/2 samples to the left (middle sample in window). voi = voi(fix(N/2)+1:length(voi);voi(length(voi)*ones(fix(N/2),1); 4、语音信号的频域分析4.1 加载“ ma1_1 ”语音数据。基于DFT 变换，画出其中一帧数据（采样频率为8kHz ，帧长为37.5ms ，每帧有300 个样点）的频域波形（对数幅度谱）。load ma1_1; x = ma1_1 (4161:4460); subplot(121), plot(x) N = 1024; k = 0:N/2-1; X = fftshift (fft (x.*hann (length (x),N); subplot(122), plot (k,20*log10 (abs(X(N/2:-1:1), axis (0 N/2-1 -inf inf ) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 10 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档4.2 编写一个用以计算信号谐波乘积谱的函数“hpspectrum.m ” ，其Matlab代码如下：function HPSx = hpspectrum(x,N,R) k = 1:R:N/2; K = length(k); X = fft (x.*hann(length(x), N); HPSx = X(k); for (r = R-1:-1:1) HPSx = HPSx.*X (1:r:r*K); end 4.3 利用函数“hpspectrum.m ” ，来估计一帧浊音的基音周期。代码如下：X = ma1_1(4161:4460); N = 1024; R = 5; HPSx = hpspectrum (x, N, R); plot (20*log10 (abs (HPSx) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 11 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档4.4 加载“vowels.mat ”语音数据， 分别画出一帧/i/ 和一帧 /u/（ 采样频率为10kHz ，帧长为30ms ，每帧有300 个样点）的基于DFT 的对数幅度谱。其Matlab 代码如下：load vowels x = vowels.i_1(2001:2300); N = 1024; k= -N/2:N/2-1; X = fftshift (fft (x.*hann (length(x),N); plot (k,20*log10(abs(X), axis(0 fix (N/2) 0 100) x = vowels.u_1(2001:2300); N= 1024; k = - N/2:N/2-1; X = fftshift (fft (x.*hann(length(x),N); plot (k,20*log10(abs(X), axis(0 fix(N/2) 0 100) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 12 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档4.5 画出一帧清音语音的基于DFT的对数幅度谱。语音数据为ma1_1中的第15701 15860 个样点（采样频率为8kHz ，帧长为20ms，每帧有160 个样点） 。load ma1_1; x = ma1-1 (4161:4460); plot (x) N = 1024; k = - N/2;N/2-1; X = fftshift (fft (x.*hann (length (x),N); plot (k,20*log10 (abs(X), axis (0 fix (N/2) - inf inf ) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 13 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档4.6加载“ timit1 ”语音数据。利用函数specgram ，画出该句语音的语谱图。其Matlab 代码如下：load timit1; NFFT = 256 ; Fs = 16000; Win = 256; Noverlap = 128; specgram(timit1, NFFT, Fs, Win, Noverlap); 4.7 加载“ gliss ”语音数据。画出该句语音的语谱图。load gliss; NFFT = 256 ; Fs = 10000; Win = 256; Noverlap = 128; specgram(gliss.i_2, NFFT, Fs, Win, Noverlap); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 14 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档4.8 基于浊音/清音分类函数voiunvoi.m和谐波乘积谱函数hpspectrum.m ，编写一个用以计算浊音短时基音周期的函数“stpitch.m ” 。其Matlab代码如下：function Fp = stpitch(x,N,Pth,Zth,NFFT,R,Fs) M = length(x); N = 2*fix(N/2); win = hann(N); Fp = zeros(M,1); voi = voiunvoi(x,N,Pth,Zth); invoi = 0; for (m=N:N/2:M) n1 = m-N+1:m-N/2; n2 = m-N/2+1:m; n = n1 n2; if (any(voi(n) if (invoi) Fp(n1) = Fmax; invoi = 0; end else HPSmax,Fmax = max(abs(hpspectrum(x(n),NFFT,R); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 15 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档if (invoi) Fp(n) = Fp(n) + win*Fmax; else Fp(n) = Fmax*ones(N/2,1) win(N/2+1:N); invoi = 1; end end end wsave = warning; warning(off); Fp = (Fp*Fs/NFFT).*(Fp./Fp); warning(wsave); 1. 调用 matlab 中的 wavread 和 wavplay 读入并播放语音文件。2. 调用 fft 命令计算语音的频谱。3. 调用 specgram 命令得到语音的语谱图。4. 调用 plot 命令绘制并观察语音的时频域波形。5. 根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时 功 率 、 短 时 过 零 率 、 基 音 频 率 的 函 数 ： function Px=stpower、 function Zx = stzerocross(x, N) 和 function HPSx = hpspectrum(x,N,R)。6. 利用短时功率和短时过零率这两个参数可以对语音信号进行浊音/清音分类。编写实现此分类功能的函数voiunvoi.m。4.9 利用函数stpitch.m ，分别求出两句语音“timit1 ”和“ timit2 ”的基音周期。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 16 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档实验二语音编码一、实验目的熟悉语音基本压缩编码的方法，观察语音压缩效果，加深对语音线性预测编码（ LPC ）的理解。二、实验内容1、编写并调试语音LPC 参数提取程序。2、编写并调试语音基音周期提取程序。3、编写并调试语音LPC 合成程序。三、实验原理语音信号中含有大量的冗余信息，采用各种信源编码技术减除语音信号的冗余度，并充分利用人耳的听觉掩蔽效应，就可以将其编码速率压缩很多倍，而仍能提供可懂语音。LPC 声码器是一种比较简单实用的语音压缩方法，其基本原理是：根据语音生成模型，将语音看作激励源通过一个线性时不变系统产生的输出，利用线性预测分析对声道参数进行估值，将求得的线性预测系数，结合基音周期等少量参数进行传输，就可以在接收端利用合成滤波器重构语音信号。线性预测系数的估计方法为：假设语音的当前样值可以用过去的M个语音样值来进行预测MiiMinxaMnxanxanxanx12121式中ia即为线性预测系数。实际值和预测值之间的均方误差可表示为nMiininxanxnE212要求均方误差总和最小，将E关于ia的偏导数设置为零，可以得到01Miininxanxknx通过采用自相关法、协方差法或格形法求解该方程，即可得到最优的ia。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 17 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档四、实验方法及程序1. 调用 xcorr 命令计算一帧语音的自相关函数。2. 调用toeplitz命令形成该帧语音的自相关矩阵。3. 调用durbin命令，采用杜宾递推算法计算该帧语音的线性预测系数。4. 编写 lpcauto.m 函数，求取一句语音信号的线性预测系数及预测残差。选择设当的窗函数对语音信号进行分幀。5. 编写 lpcpitch 函数，由残差信号计算该句语音的基音周期。6. 编写 lpcgain 函数，由预测残差能量，求出该句语音的增益。7. 编写 lpcsyn 函数， 由该句语音的基音周期、预测残差能量和增益进行LPC 合成。五、实验结果与分析1. 如 何对 全极点模型的线性预测参数进行提取？加窗等预处理对提取结果有何影响？自相关法、协方差法求解线性预测方程组、格型法。协方差法是先把计算均方误差的间隔N固定下来，再研究它对计算n（ i,1 ）的其中N是基音周期值。添加汉明窗和汉宁窗可以得到效果不错的波形，很明显两个波形相似，矩形窗的波形相对不是很平稳。通过改变帧长、帧移后波形更加的清晰，特征更加的明显，更方便观察。2. 如何实现线性预测参数的各种表现方式之间的转换？对于语音信号，确定了各线性预测系数后，根据H（ z）可得其频率响应的估值即 LPC 谱。3. 使用预测残差求基音周期有何优点？该方法能进一步去除声道共振峰的影响，简化清浊音判决过程，提高判决精度，有效降低发生基音倍频和半频的错误。4. 如何根据线性预测系数求得LPC 频谱？对于语音信号，确定了各线性预测系数后，根据H（ z）可得其频率响应的估值即 LPC 谱。5. LPC 频谱与基于DFT 求得的语音对数幅度谱相比有何不同？LPC 谱优点在于可以很好地表示共振峰结构而不出现额外的蜂起和起伏，基本不含有噪声，其频率灵敏度与人耳不相匹配；DET 得到的对数幅度谱受基频谐波的影响，最大值只出现在谐波频率上，共振峰测定误差较大。6. LPC 合成语音的质量如何？有何改进措施？LPC 合成技术的优点是简单直观。其合成过程实质上只是一种简单的解码和拼接过程。另外，由于波形拼接技术的合成基元是语音的波形数据，保存了语音的全部信息，因而对于单个合成基元来说能够获得很高的自然度。但是，由于自然语流中的语音和孤立状况下的语音有着极大的区别，如果只是简单地把各个孤立的语音生硬地拼接在一起，其整个语流的质量势必是不太理想的。而LPC 技术从本质上来说只是一种录音+重放，对于合成整个连续语流LPC 合成技术的效果是不理想的。因此， LPC 合成技术必须和其他技术相结合，才能明显改善LPC 合成的质量。实验过程记录与结果分析1、计算短时自相关函数1.1分别计算正弦信号和白噪声的短时自相关函数，估计正弦信号的基音周期。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 18 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档其 Matlab 代码如下：x=sin(2*pi*0.01*(0:499); r,eta=xcorr(x,100,unbiased); stem(eta,r); w=randn(500,1); r,eta=xcorr(w,100,unbiased); stem(eta,r); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 19 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档1.2计算正弦信号叠加白噪声的短时自相关函数，试估计正弦信号的基音周期。其 Matlab 代码如下：x=sin(2*pi*0.01*(0:499); w=randn(500,1); x1=x+w; r,eta=xcorr(x1,100,unbiased); stem(eta,r); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 20 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档1.3 分别画出一帧浊音和一帧清音的语音自相关波形（采样频率为10kHz ，帧长为 25.6ms ，每帧有256 个样点） ，估计浊音的基音周期。load digits; x=digits.three1; m=2756; N=256; n=m-N+1:m; r,eta=xcorr(x(n),250,unbiased); plot(eta,r); m=500; N=256; n=m-N+1:m; r,eta=xcorr(x(n),250,unbiased); plot(eta,r); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 21 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档2、求取 LPC系数2.1加载“ digits.three1”语音数据。基于自相关法，求出其中一帧数据（采样频率为10kHz ，帧长为25.6ms ，每帧有256 个样点）的14 阶 LPC 系数。画出相应的 LPC 谱。load digits; x=digits.three1; m=2756; N=256; n=m-N+1:m; M=14; r,eta=xcorr(x(n),250,unbiased); Rx=toeplitz(r(M+1:2*M); rx=r(M+2:2*M+1); a=Rxrx; NFFT=1024; k=1:NFFT/2; X=fft(x(n).*hann(N),NFFT); Theta=1./fft(1;-a,NFFT); plot(k,20*log10(abs(353*Theta(k) X(k) axis(0 NFFT/2 inf inf) 2.2 编写一个用以实现杜宾递推算法的函数“durbin.m ” ，其 Matlab 代码如下：function a,xi,kappa=durbin(r,M) kappa = zeros (M,1); a = zeros (M, 1); xi=r(1); zeros(M,1); for (j=1:M) kappa(j)=r(j+1)-a(1:j-1)*r(j:-1:2)/xi(j); a(j)=kappa(j); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 22 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档a(1:j-1)=a(1:j-1)-kappa(j)*a(j-1:-1:1); xi(j+1)=a(1:j-1)-kappa(j)*a(j-1:-1:1); xi(j+1)=xi(j)*(1-kappa(j)2); end 2.3 利用函数“durbin.m ” ，计算2.1 中语音数据帧的14 阶 LPC 系数，并与2.1中的结果进行比较：x=digits.three1; m=2756; N=256; n=m-N+1:m; M=14; r,eta=xcorr(x(n),250,unbiased); aLD,xi,kappa=durbin(r(M+1:2*M+1),M); a, aLD, norm(a-aLD) a = 0.3512 0.2072 0.3808 0.3163 0.0522 -0.1639 -0.2903 -0.2898 -0.1222 0.1960 0.3492 0.3024 0.1713 -0.8258 aLD = 0.3512 0.2072 0.3808 0.3163 0.0522 -0.1639 -0.2903 -0.2898 -0.1222 0.1960 0.3492 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 23 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档0.3024 0.1713 -0.8258 ans = 1.0047e-14 2.4 编写一个用以实现反射系数转换为LPC 系数的函数“rf2lpc.m ” ，其Matlab代码如下：function a = rf2lpc(kappa) M = length(kappa); a = zeros(M,1); for (j=1:M) a(j) = kappa(j); a(1:j-1) = a(1:j-1) - kappa(j)*a(j-1:-1:1); end 2.5 编写一个用以实现LPC 系数转换为反射系数的函数“lpc2rf.m ” ，其 Matlab代码如下：function kappa = lpc2rf(a) M = length(a); kappa = zeros(M,1); for (j=M:-1:1) kappa(j) = a(j); a(1:j-1) = (a(1:j-1) + a(j)*a(j-1:-1:1)/(1 - kappa(j)2); end 2.6 使用函数“rf2lpc.m ”和“ lpc2rf.m ” ，检验反射系数和LPC 系数相互转换结果。norm(kappa-lpc2rf(aLD) norm(aLD-rf2lpc(kappa) ans = 1.4760e-14 ans = 0 3、语音信号的逆滤波3.1 利用2.3 中求出的一帧语音数据的LPC 系数，构造逆滤波器，并画出该帧语音信号的残差波形。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 24 页，共 30 页 - - - - - - - - - - 资料收集于网络如有侵权请联系网站删除谢谢精品文档ehat=filter(1;-1,1,x(n); plot(x(n) ehat) 3.2 编写一个用以实现语音信号LPC 分析的函数“lpcauto ” ，其Matlab代码如下：function ar,xi,e,m = lpcauto(x,M,win,Olap) Nx = length(x); N = length(win); if (N = 1) N = win; win = ones(N,1); end F = fix(Nx-Olap)/(N-Olap); ar = zeros(M+1,F); xi = zeros(M+1,F); e = zeros(Nx,1); m = zeros(F,1); n = 1:N; n1 = 1:Olap; n2 = N-Olap+1:N; n3 = Olap+1:N; win1 = win(n1).