语音信号的短时自相关序列求解以及xcorr与autocorr问题

对于语音来说，短时自相关函数分析是一个重要的方法，能够用来求得浊音的基音周期，也可以用来求得语音识别中的特征参数。它的短时自相关函数为：

         

但是，在计算短时自相关时，窗选语音段为有限长度N，而求和上限为N-1-k，因此当k增加时可用于计算的数据就越来越少了，从而导致k增加时自相关函数的幅度减小。为了解决这个问题，提出了语音修正的短时自相关。修正的短时自相关函数，其定义如下：

     为了用Matlab计算一帧语音的自相关序列，可以采用Matlab中提供的xcorr和autocorr函数，也可以自己按照上面的公式动手编写。

代码如下：

加上frame1帧长为256的信号。

%用修正的自相关函数计算自相关，采用两个不同长度的窗口

N=128; %
K=128;
A=[];
for k=1:K
    sum=0;
    for m=1:N
        sum=sum+frame1(m)*frame1(m+k-1);
    end
    A(k)=sum;
end
for k=1:K
    A1(k)=A(k)/A(1);
end
subplot(4,1,2);
plot(A1);
title('修改的自相关函数求自相关');

该方法计算得出的自相关序列只有128个点，但是如果想得到更多的点呢？ 可是用Matlab提供的函数xcorr和autocorr可以计算数量为帧长256个自相关序列，其中xcorr计算出来的为对称序列,需要去掉负延迟(lag)的自相关序列。而autocorr只取了一边，不是对称的。

下面介绍下xcorr与autocorr的具体区别：

1.xcorr

C=xcorr(A,B),求向量A与B的互相关系数。如果A和B都是长度为M的向量，则返回2*M-1个胡相关系数C。
如果A和B的长度不一样，短的向量补0，然后计算互相关。如果A是行向量，C也是行向量。

xcorr(A),如果A是一个向量，则求自相关序列。如果A是一个M*N的矩阵，则求出的结果为2M-1行N^2列的自相关序列。延迟为0
的点位于该序列的中间。


xcorr(A,MaxLag),MaxLag=M-1,M为向量A的长度。求延迟为-MaxLag到MaxLag之间的相关函数。


[c,lags]=xcorr(A,MaxLag,ScaleOpt);返回的lags为延迟下标


ScaleOpt为归一化选项：'biased':自相关序列乘以1/M ;‘unbiased’：自相关序列乘以 1/(M-abs(lags))
 'coeff':归一化序列让延迟为0的自相关序列为1.
‘none’：不归一化


例子： s=[1  2  3];  
r=xcorr(s);
 r =3.0000    8.0000   14.0000    8.0000    3.0000
上面的矩阵，M=3,最后得到5个结果，其中第三个是自己和自己相乘，最后相加的结果，值最大1*1+2*2+3*3=14。而第二个和
第四个分别是间隔正负1的结果也就是1*2+2*3=8，2*1+3*2=8。第1个和第五个分别是间隔正负2，也就是1*3=3，3*1=3。

可见：xcorr默认的MaxLag=M-1=2.ScaleOpt默认为none.
MaxLag=2;
r=xcorr(s,MaxLag,'biased');
r =1.0000    2.6667    4.6667    2.6667    1.0000
相当于乘以1/3


MaxLag=2;
[r,lags]=xcorr(s,MaxLag,'unbiased');
 r=3.0000    4.0000    4.6667    4.0000    3.0000
abs(lags)= 2 1 0 1 2
相当于乘以1/(M-abs(lags))，第一个为1/(3-2)=1; 1/(3-1); 1/(3-0);1/(3-1);1/(3-2)

MaxLag=2;
[r,lags]=xcorr(s,MaxLag,'coeff'); (注该方法相当于修正的自相关函数)
 r=0.2143    0.5714    1.0000    0.5714    0.2143
相当于r=r./r(M)  即3/14=0.2143    8/14=0.5714,  14/14=1.00;

2.autocorr
[acf,lags,bounds] = autocorr(y) 
[acf,lags,bounds] = autocorr(y,numLags,numMA,numSTD)
计算随机时间变量y的自相关函数，当没有输出参数，将画出在置信区间的自相关函数。
 numLags：为自相关函数的延迟，默认为min[20,length(y)-1];该函数忽略了延迟小于0的序列。

下面用修正的自相关函数和xcorr,autocorr来计算一帧语音的自相关序列

代码如下：

[plain]  view plain copy

1. clc;  
2. clear all;  
3. [y,fs,nbits]=wavread('1b.wav');  
4. frameSize=256;  
5. overLap=128;  
6. lag=255;  
7. frameMat=buffer2(y,frameSize,overLap);  
8. [ndim,nFrame]=size(frameMat);  
9. frame1=frameMat(:,50);  
10. subplot(4,1,1);  
11. %plot([1:frameSize]/fs,frame1);  
12. %画出该语音帧的时域波形图  
13. plot([1:frameSize],frame1);  
14.   
15. %用修正的自相关函数计算自相关，采用两个不同长度的窗口  
16. N=128;  
17. K=128;  
18. A=[];  
19. for k=1:K  
20.     sum=0;  
21.     for m=1:N  
22.         sum=sum+frame1(m)*frame1(m+k-1);  
23.     end  
24.     A(k)=sum;  
25. end  
26. for k=1:K  
27.     A1(k)=A(k)/A(1);  
28. end  
29. subplot(4,1,2);  
30. plot(A1);  
31. title('修改的自相关函数求自相关');  
32. %xcorr求自相关时，没有减去均值，没有归一化，而autocorr求自相关是减去了均值也做了归一化处理  
33. %%%%%xcorr不减去均值求自相关  
34.   
35. [c,lags]=xcorr(frame1,lag, 'unbiased'); %unbiased估计不会衰减,默认的最大lag为frameSize-1;  
36. %d=c./c(lag+1);%归一化，除以能量 已经用1/(M-abs(lags))进行了归一化  
37. subplot(4,1,3);  
38. plot(lags(lag+1:end),c(lag+1:end));  
39. %axis([0,256,-1,1]);  
40. title('xcorr求自相关unbiased(没有减去均值)');%该函数实现本身没有减掉均值做相关  
41.   
42.   
43.   
44.   
45. [c,lags]=xcorr(frame1,lag, 'biased'); %unbiased估计不会衰减,默认的最大lag为frameSize-1;  
46. %d=c./c(lag+1);%归一化，除以能量 已经用1/(M-abs(lags))进行了归一化  
47. subplot(4,1,4);  
48. plot(lags(lag+1:end),c(lag+1:end));  
49. %axis([0,256,-1,1]);  
50. title('xcorr求自相关biased(没有减去均值)');%该函数实现本身没有减掉均值做相关  
51.   
52. [c,lags]=xcorr(frame1,lag, 'coeff');   
53. %d=c./c(lag+1);%归一化，除以能量 已经用1/(c(lag+1))进行了归一化  
54. figure;  
55. subplot(4,1,1);  
56. plot(lags(lag+1:end),c(lag+1:end));  
57. %axis([0,256,-1,1]);  
58. title('xcorr求自相关coeff(没有减去均值)');%该函数实现本身没有减掉均值做相关  
59.   
60. %%%%%xcorr减去均值求自相关  
61. frame2=frame1-mean(frame1);%减掉均值  
62. [c1,lags1]=xcorr(frame2,lag, 'unbiased'); %unbiased估计不会衰减,默认的最大lag为frameSize-1;  
63. %d1=c1./c1(lag+1);%归一化，除以能量  
64. %figure;  
65. subplot(4,1,2);  
66. plot(lags1(lag+1:end),c1(lag+1:end));  
67. %axis([0,256,-1,1]);  
68. title('xcorr求自相关unbiased(减去均值)');%该函数实现本身没有减掉均值做相关  
69.   
70.   
71.   
72.   
73. [c1,lags1]=xcorr(frame2,lag, 'coeff'); %unbiased估计不会衰减,默认的最大lag为frameSize-1;  
74. %d1=c1./c1(lag+1);%归一化，除以能量  
75. %figure;  
76. subplot(4,1,3);  
77. plot(lags1(lag+1:end),c1(lag+1:end));  
78. %axis([0,256,-1,1]);  
79. title('xcorr求自相关coeff(减去均值)');%该函数实现本身没有减掉均值做相关  
80.   
81.   
82. %%%%用auotcorr求自相关  
83. [c2,lags2,bound]=autocorr(frame1,lag);  
84. subplot(4,1,4);  
85. plot(lags2(1:end),c2(1:end));  
86. %axis([0,256,-1,1]);  
87. title('autocorr求自相关');%该函数减掉了均值在做相关,最后做了归一化  

如下图所示：

从上图可以看出autocorr(y)求得的结果相当于xcorr(y-mean(y),MaxLag,'coeff')的结果。