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

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

南京邮电大学通达学院

实 验 报 告

实验名称 语音时频参数的提取和应用实验

语音编码实验

课程名称 信息处理技术专业综合实验

班级学号 姓 名

开课时间 2016/2017 学年, 第 二 学期

实验一 语音时频参数的提取和应用

一、实验目的

熟悉基本的语音时频参数提取方法，观察并比较各类参数采用不同提取方法的效果，了解这些时频参数在语音处理中的应用。

二、实验内容

1、编写并调试语音时域参数提取程序，包括短时功率、短时过零率等函数。

2、编写并调试语音频域参数提取程序，包括频谱、语谱图、基音频率、共振峰等。 3、利用提取出的时频参数对语音信号进行静/清/浊音的划分。

三、实验原理

语音信号具有短时平稳性，需对其进行分段（帧）处理（帧长一般取10至30ms之间）。 根据激励方式的不同，语音信号可以分成浊音和清音两大类：

浊音：激励源可模拟成准周期性的脉冲串。该周期称为基音周期，其倒数就是基音频率，基音频率一般在几百赫兹以内。 清音：激励源可模拟成随机白噪声。 简化的语音生成数学模型如图所示：

基音周期准周期脉冲序列发生器浊音时变线性系 统清音随机噪声发生器语音取样值

对语音时\\频波形进行绘制与观察，有助于了解语音的清/浊音等特性。对基音周期等声门参数和共振峰等声道参数进行提取与分析，是语音处理与应用的重要内容。

对于时变语音信号x(n)，可以计算它的短时功率P(m)和短时过零率Zxx增益控制时变参数(m)。第m帧（帧长为N

个样点）的功率计算公式为：

1Px(m)NnmN1x(n)

m2过零率计算公式为：

1Zx(m)NnmN1msign{x(n)}sign{x(n1)}2

对于浊音语音，可以利用其频谱X()具有丰富的谐波分量的特点，求出其谐波乘积谱，计算公式为：

HPSx()X(r)r1R

式中，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，第一帧0~100Hz有10个峰值点，基频约为10Hz，第二帧0~100Hz有5个峰值点，基频约为20Hz，第三帧0~100Hz有5个峰值点，基频约为20Hz， 第四帧0~100Hz有4个峰值点，基频约为25Hz。4. 如何根据基于DFT的对数幅度谱估计某幀浊音的基音周期？ 0~1172为一个基因周期。

5. 如何根据基于DFT的对数幅度谱，估计出共振峰频率？

有7个峰点，因此基频约为167。

6. 时域对语音信号进行加窗，反映在频域，其窗谱对基于DFT的对数幅度谱有何影响？如何估计出窗谱的主瓣宽度？

因为时间窗幅度的傅立叶变换为脉冲状态,则频谱偏差很小,当定常过程x(n)的平均值β为零且N很大时 E[TN(ω)]=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 % 十个孤立字的发音 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

447288 double

2.3利用命令“length”和“音信号timit1的长度和维数：m = length(timit1) [m,n] = size(timit1) （粘贴实验结果） m =

55911 m =

55911

55911x1 size”查看语

n = 1

2.4通过以下命令，分别查看语音信号timit1中间的1个、2个和20个数据值：

timit1(5001) timit1([5001,5003]) timit1(5001:5020) （粘贴实验结果） ans =

-0.06 ans =

-0.06 0.0111

ans =

-0.06 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 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.

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;

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))

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 = (Px>Pth*max(Px)) & (Zxvoi s(fix(N/2),1)];

4、语音信号的频域分析

=

[voi(fix(N/2)+1:length(voi));voi(length(voi))*one

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), (abs(X(N/2:-1:1)))),

axis ([0 N/2-1 -inf inf ])

plot

(k,20*log10

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)))

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])

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 ])

4.6 加载“timit1”语音数据。利用函数specgram，画出该句语音的语谱图。其Matlab代码如下：

load timit1;

NFFT = 256 ; Fs = 16000; Win = 256; Noverlap = 128;

specgram(timit1, Noverlap);

NFFT, Fs, Win,

4.7加载“gliss”语音数据。画出该句语音的语谱图。

load gliss;

NFFT = 256 ; Fs = 10000; Win = 256; Noverlap = 128;

specgram(gliss.i_2, NFFT, Fs, Win, Noverlap);

4.8 “stpitch.m”。其Matlab代码如下： function M = length(x); N = 2*fix(N/2); win = hann(N); Fp = zeros(M,1);

voi = voiunvoi(x,N,Pth,Zth); invoi = 0;

Fp

=

stpitch(x,N,Pth,Zth,NFFT,R,Fs)

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)));

if (invoi) Fp(n) = Fp(n) + win*Fmax; else Fp(n) win(N/2+1:N)];

invoi = 1; end end

=

Fmax*[ones(N/2,1)

=

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”的基音周期。

实验二 语音编码

一、实验目的

熟悉语音基本压缩编码的方法，观察语音压缩效果，加深对语音线性预测编码（LPC）的理解。

二、实验内容

1、编写并调试语音LPC参数提取程序。 2、编写并调试语音基音周期提取程序。 3、编写并调试语音LPC合成程序。

三、实验原理

语音信号中含有大量的冗余信息，采用各种信源编码技术减除语音信号的冗余度，并充分利用人耳的听觉掩蔽效应，就可以将其编码速率压缩很多倍，而仍能提供可懂语音。 LPC声码器是一种比较简单实用的语音压缩方法，其基本原理是：根据语音生成模型，将语音看作激励源通过一个线性时不变系统产生的输出，利用线性预测分析对声道参数进行估值，将求得的线性预测系数，结合基音周

期等少量参数进行传输，就可以在接收端利用合成滤波器重构语音信号。 线性预测系数的估计方法为：

假设语音的当前样值可以用过去的M个语音样值来进行预测

~xna1xn1a2xn2aMxnMaixnii1M

式中a即为线性预测系数。实际值和预测值

i之间的均方误差可表示为

要求均方误差总和最小，将E关于a的偏导

iE2nxnaixninni1M2数设置为零，可以得到

Mxnkxnaxni0ini1

通过采用自相关法、协方差法或格形法求解

该方程，即可得到最优的a。

i

四、实验方法及程序

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 分别计算正弦信号和白噪声的短时自相关函数，估计正弦信号的基音周期。其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);

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);

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);

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=Rx\\rx;

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);

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.28 -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.28 -0.1222 0.1960 0.3492 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系数，构造逆滤波器，并画出该帧语音信号的残差波形。

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)./(win(n1)+win(n2)+eps); win2 = win(n2)./(win(n1)+win(n2)+eps);

for (f=1:F)

[r,eta] = xcorr(x(n).*win,M,'biased'); [a,xi(:,f),kappa] durbin(r(M+1:2*M+1),M);

ar(:,f) = [1; -a];

ehat = filter(ar(:,f),1,x(n)); e(n)

m(f) = n(N); n = n + (N-Olap); end

3.3加载“timit1”语音数据，利用“lpcauto”函数对该句语音进行LPC分析。画出语音信号及其残差波形。

load timit1; x=timit1; M=14; N=256;

[ar,xi,e,m]=lpcauto(x,M,hann(N),N/2); subplot(121), plot(x) subplot(122),plot(e) axis([0,60000,-1,1])

=

[e(n(n1)).*win2

+

ehat(n1).*win1; ehat(n3)]; % Overlap-add.

=

plot([x,e]) soundsc(x) soundsc(e) soundsc(x-e)

4、LPC谱估计

4.1编写一个用以计算和显示LPC谱的函数“lpcplot”，其Matlab代码如下：

function lpcplot(A,Nfft,Fs,m) [M,N] = size(A); if (N==1)

[Theta,F] = freqz(1,A,Nfft,Fs); plot(F,20*log10(abs(Theta)));

xlabel('Frequency, {\\it F} [Hz]'); ylabel('Magnitude, [dB]');

else

if (length(m) ~= N)

error('The column dimension of A must be equal to the length of m.')

end

Theta = zeros(Nfft,N); for (n=1:N) [Theta(:,n),F] freqz(1,A(:,n),Nfft,Fs);

end

MeshHndl

meshz(m,F,20*log10(abs(Theta)));

axis

ij;

view(-45,45);

set(MeshHndl,'MeshStyle','Column');

axis tight; axis 'auto y'; axis 'auto z';

= =

|\heta(\\omega)|

xlabel('Sample Number, {\\it n}'); ylabel('Frequency, {\\it F} [Hz]'); zlabel('Magnitude, [dB]');

end

4.2利用“lpcplot”函数，画出3.3中语音数据的几帧（采样频率为16kHz，帧长为32ms，每帧有512个样点）LPC谱。

plot(x(m(71):m(147)))

|\heta(\\omega)|

lpcplot(ar(:,71:147),512,16000,m(71:147))

lpcplot(ar(:,82),512,16000,0)

lpcplot(ar(:,110),512,16000,0)

lpcplot(ar(:,134),512,16000,0)