Matlab用xcorr求自相关函数出现的一点问题

 用Matlab做随机信号分析实验仿真的时候出现一个问题：使用xcorr计算信号的自相关函数并求Rx(0)=max(Rx)，然后通过E(x.^2)求均方值    σx^2，理论二者都是表示平均功率，应该相等，但是实验过程中得出高斯白噪声的均方值0.1495,而Rx(0)可达30以上。如下图：

经过阅读matlab的help文档，发现原因是xcorr缺省计算原始的没有标准化的相关性（求和而非求均值），需要加上’unbiased’选项计算自相关函数的无偏估计。

==========================================================================================================================

以下是整理的xcorr函数原始help文档和部分翻译，供参考：

xcorr

Cross-correlation

Syntax

c=xcorr(x,y)
c=xcorr(x)
c=xcorr(x,y,'option')
c=xcorr(x,'option')
c=xcorr(x,y,maxlags)
c=xcorr(x,maxlags)
c=xcorr(x,y,maxlags,'option')
c=xcorr(x,maxlags,'option')
[c,lags]=xcorr(...)

Description

xcorr estimates the cross-correlation sequence of a random process.Autocorrelation is handled as a special case.

Xcorr估计一个随机过程的互相关序列，自相关是当做一个特例处理的。

The true cross-correlation sequence is

真实的互相关序列是

                    

where x_n and y_n are jointly stationaryrandom processes, −∞ <n<∞, and E {·} is the expected value operator. xcorr mustestimate the sequence because, in practice, only a finite segment of onerealization of the infinite-length random process is available.

当xn和yn是....随机过程，−∞ < n <∞, E {·}是期望值算子，xcorr必须估计序列，因为在实践中，无限长随机过程的一个实现只有一个有限段是可用的。

c=xcorr(x,y) returns the cross-correlation sequence in alength 2*N-1 vector, where x and y are lengthN vectors (N>1).If x and y are not the same length, the shorter vector is zero-padded to thelength of the longer vector.

By default, xcorr computes raw correlations with no normalization.

缺省情况下xcorr计算原始的没有标准化的相关性（这里就是开始错误的原因所在，请注意以下表达式，貌似下面计算的是确定性信号的自相关函数，不同于随机信号，求指正）

                    

The output vector c has elements given by c(m) = R_xy(m−N),m=1, ..., 2N−1.

输出向量c的元素通过c(m) = R_xy(m−N),m=1, ..., 2N−1.给出

In general, the correlation function requires normalizationto produce an accurate estimate (see below).

一般情况下，相关性函数需要标准化以得出一个精确的估计（见下面）

c = xcorr(x) is the autocorrelationsequence for the vector x. If x is anN-by-P matrix, c is amatrix with 2N-1 rows whoseP² columns contain thecross-correlation sequences for all combinations of the columns of x. For moreinformation on matrix processing with xcorr, seeMultiple Channels.

c = xcorr(x)得出向量x的自相关序列，如果x是N*P的矩阵，c是一个2N-1行的矩阵，它的P^2列包含x的所有列组合的互相关序列。xcorr矩阵处理的更多信息，请参阅多Multiple Channels.

xcorr produces correlations identically equal to 1.0 at zero lag only whenyou perform an autocorrelation and only when you set the 'coeff' option. Forexample,

只有当你求一个自相关和只有当你设置了'coeff“选项时，xcorr在零滞后产生恒等于1.0的相关性。例如，

x=0:0.01:10;X = sin(x);[r,lags]=xcorr(X,'coeff');  max(r)

c = xcorr(x,y,'option') specifies anormalization option for the cross-correlation, where'option' is

c = xcorr(x,y,'option')为互相关指定一个标准化的选项，当选项是

• 'biased': Biased estimate of the cross-correlation function
• 'biased':互相关函数的有偏估计（这里涉及有偏和无偏估计，有兴趣可以参考相关文献，此处不仔细讲解）

                        

• 'unbiased': Unbiased estimate of the cross-correlation function
• 'unbiased':互相关函数的无偏估计

                     

• 'coeff': Normalizes the sequence so the autocorrelations at zero lag are identically 1.0.
• 'coeff':标准化序列使自相关在零滞后等于1.0
• 'none', to use the raw, unscaled cross-correlations (default)
• 'none'，使用原始的未缩放的互相关性（缺省）。

See [1] for moreinformation on the properties of biased and unbiased correlation estimates.

c = xcorr(x,'option') specifies one of theabove normalization options for the autocorrelation.

c = xcorr(x,y,maxlags) returns thecross-correlation sequence over the lag range [-maxlags:maxlags]. Output c haslength 2*maxlags+1.

c = xcorr(x,maxlags) returns theautocorrelation sequence over the lag range [-maxlags:maxlags]. Output c haslength 2*maxlags+1. If x is anN-by-P matrix, c is a matrix with2*maxlags+1 rows whoseP² columns contain the autocorrelationsequences for all combinations of the columns of x.

c = xcorr(x,y,maxlags,'option')specifies both a maximum number of lags and a scaling option for thecross-correlation.

c = xcorr(x,maxlags,'option') specifiesboth a maximum number of lags and a scaling option for the autocorrelation.

[c,lags] = xcorr(...) returns a vector of thelag indices at which c was estimated, with the range [-maxlags:maxlags]. Whenmaxlags is not specified, the range of lags is [-N+1:N-1].

In all cases, the cross-correlation or autocorrelation computed by xcorrhas the zeroth lag in the middle of the sequence, at element or row maxlags+1(element or row N if maxlags is not specified).

Examples

The second output, lags, is useful for plotting the cross-correlation orautocorrelation. For example, the estimated autocorrelation of zero-meanGaussian white noisec_ww(m) can be displayed for -10 ≤m ≤ 10 using:

ww = randn(1000,1);[c_ww,lags] = xcorr(ww,10,'coeff');stem(lags,c_ww)

Swapping the x and y input arguments reverses (and conjugates) the outputcorrelation sequence. For row vectors, the resulting sequences are reversedleft to right; for column vectors, up and down. The following example illustratesthis property (mat2str is used for a compact display of complex numbers):

x = [1,2i,3]; y = [4,5,6];[c1,lags] = xcorr(x,y);c1 = mat2str(c1,2), lagsc2 = conj(fliplr(xcorr(y,x)));c2 = mat2str(c2,2)

For the case where input argument x is a matrix, the output columns arearranged so that extracting a row and rearranging it into a square arrayproduces the cross-correlation matrix corresponding to the lag of the chosenrow. For example, the cross-correlation at zero lag can be retrieved by:

randn('state',0)X = randn(2,2);[M,P] = size(X);c = xcorr(X);c0 = zeros(P); c0(:) = c(M,:)    % Extract zero-lag row

You can calculate the matrix of correlation coefficients that the MATLABfunctioncorrcoef generates by substituting:

c = xcov(X,'coef')

in the last example. The function xcov subtractsthe mean and then calls xcorr.

Use fftshift to move the second half of the sequence starting at the zeroth lag to thefront of the sequence. fftshift swaps the first and second halves of a sequence.

Algorithm

For more information on estimating covariance and correlation functions,see[1].

References

[1] Orfanidis, S.J., Optimum Signal Processing. AnIntroduction. 2nd Edition, Prentice-Hall, Englewood Cliffs, NJ, 1996.