目录
[隐藏]
1 定义1.1 统计学
1.2 信号处理
2 自相关函
3 自相关函数举例
4 参考文献
定义[编辑]
自相关函数在不同的领域,定义不完全等效。在某些领,自相关函数等于自协方
统计学[编辑]
将一个有序的随机变量系列与自身相比较,这就是自相关函数在统计学中的定义。每不存在相位差的系列,都与其自相似,在此况下,自相关函值最大。如果系列中的组成部分相之间存在相关性(不再是随机的),则由以下关值方程计的值不再为零,这样的组成部分
......... 期望值。
........ 在t(i)时的随机变量
........ 在t(i)时的预期
.... 在t(i+k)的随机变量
.... 在t(i+k)时的预期值。
......... 为方差。
所得的自相关值R的取值范围为[-1,1],1为最大正相关值,-1则为最大负相关,0为不
信号处理[编辑]
在信息分析中,通常将自相关函数称之为自协方差方。 用来描述信在不同
的,信息函数值
,其中“
”是卷积算符,
为取共轭。
自相关函数的性
以下以一维自相关函数为例说明其性质,多维的情况方便地从一维情推广得
对称性:从定义显然可以看出R(i)?=?R(?i)。连续型自相函数为偶
当f为实函数
当f是复函数时,该自相关函数厄米函数,满
其中星号表示共轭。
连续型实自相关函数的峰值在原点取得,即于任何延时 τ,
。该结论可直接有柯西-施瓦茨不等式得到。离型自相关函数有此结
周期函数的自相关函数是具有与原数相同周期的函
两个相互无关的函数(即对于所有 τ,两函的互相关均为0)之和自相关函数等于各自相关函数
由于自相关函数是一种特殊的互相关函数,所它具有后者所有性
连续时间白噪声信号的自相关函数是个δ函数,在除 τ?=?0 之外的
维纳-辛钦定理(Wiener–Khinchin theorem)表明,自相函数和功率谱密度函数是对傅里叶变
实值、对称的自相关函数具有实称的变换函,因此此时维纳-辛钦定中的复指数项可以写成下的余弦
自相关函数举
白噪声的自相关函
参考文献[编辑]
自相关函数
自相关函数在不同的领域,定义不完等效。在某些域,自相关函数等同于自协
统计学
R(k) = \frac{E[(X_i - \mu)(X_{i+k} - \mu)]}{\sigma^2}
信号处理
R_f(\tau) = f(\tau) * f^*(-\tau)= \int_{-\infty}^{\infty}
f(t+\tau)f^*(t)\, dt = \int_{-\infty}^{\infty} f(t)f^*(t-\tau)\,
dt,其中“*”是卷积算符,(\cdot)^*为取
同一时间函数瞬时t和t+a的两个值相乘积平均值作为延迟时间t函数,它是信号与延迟后信号之间相似性的量。延迟时间为零时,则成为信号的
编辑本段
自相关函数的性质
以下以一维自相关函数为例说明性质,多维的情况可便地从一维情况
对称性:从定义显然可以看出R(i) = R(?i)。连续型自相关函
当f为实函数
R_f(-\tau) = R_f(\tau)\,
当f是复函数时,该自相关函数是
R_f(-\tau) = R_f^*(\tau)\,
其中星号表示
连续型实自相关数的峰值在原点取得,即对于任何时 τ,均有 |R_f(\tau)| \leq R_f(0)。该论可直接有柯西-施瓦兹不等式得到。散自相关函数亦有此结
周期函数的自相关函数是具有与原函
两个相互无关的函数(即对于有 τ,两数的互相关均为0)之和自相关函数等于各自自
由于自相关函数是一种特殊的互相关函数,所它具有后者的所
连续时间白噪声信号的自相关函数是个δ函数,在除 τ = 0 之外的所有
维纳-辛定理(Wiener–Khinchin theorem)表明,自相函数和功率谱密度函数是一
R(\tau) = \int_{-\infty}^\infty S(f) e^{j 2 \pi f \tau} \, df
S(f) = \int_{-\infty}^\infty R(\tau) e^{- j 2 \pi f \tau} \, d\tau.
实值、对称的自相关函数具有实对称的变函数,因此此维纳-
定理中的复指数项可以写成如
R(\tau) = \int_{-\infty}^\infty S(f) \cos(2 \pi f \tau) \, df
S(f) = \int_{-\infty}^\infty R(\tau) \cos(2 \pi f \tau) \, d\tau. 编辑
自相关函数举例
白噪声的自相关函数
r_{nn} = \mathbb{E} \{ n(t) n(t-\tau) \} = \delta ( \tau )
FFT实现自相关函数
FFT 实现
N=38;
noise=(randn(1,N)+1i*randn(1,N))/sqrt(2); f1=0.11;
f2=0.15;
f3=0.23;
SNR1=20;
SNR2=18;
SNR3=17;
A1=10^(SNR1/20);
A2=10^(SNR2/20);
A3=10^(SNR3/20);
signal1=A1*exp(1i*2*pi*f1*(0:N-1));
signal2=A2*exp(1i*2*pi*f2*(0:N-1));
signal3=A3*exp(1i*2*pi*f3*(0:N-1));
un=signal1+signal2+signal3+noise;
Uk=fft(un,2*N);
Sk=(1/N)*abs(Uk).^2;r0=ifft(Sk);
r1=[r0(N+2:2*N),r0(1:N)];
r=xcorr(un,N-1,'biased');
r11=real(r1);
r12=imag(r1);
r1=real(r);
r2=imag(r);
m=1-N:N-1;
subplot(2,2,1);
stem(m,r11,'o');
xlabel('m');
ylabel('
title('基于 FFT 的
subplot(2,2,2);
stem(m,r12,'o');
xlabel('m');
ylabel('
subplot(2,2,3);
stem(m,r1);
xlabel('m');
ylabel('
title('基于直接计算的
subplot(2,2,4);
stem(m,r2);
xlabel('m');
ylabel('虚部
');
自相关函数
信号检测与估计
作业要求
1、产生两个点数为2000的高斯分布随机数列,并画出其二维概率密度; 2、用给定的两个录音,以采频率22050进行采,并分别画出其自相关函数曲线; 3、自己任选一个信号,最好是自己录制的信号,采频率自己确定,画出其自相关
一、画二维概
首先需要产生两个点数为2000的高斯分布随机数,以用matlab自带的函数randn()。后为了计出二维概率密度,我们需要知道五参数,两个高斯数列的均值、方差以相系数。可以用normfit命令计算均吗,然后通过cov(x,y)求出协方矩阵。然后就可以通过mvnpdf命令计算出对应坐标点的概率密函数值。为我们是要画一个二维图,所以需要事先打好网格,我将网格的范围定为[-5,5]*[-5,5],步为0.1。然后计算每一个网格点的概率度函数值,最后即
程序代码如下:
clear
x=randn(2000,1);%产生2000个标准正态布随机
y=randn(2000,1);%产生2000个标准正态布随机
[mulx,sigmax]=normfit(x);%计y1均
[muly,sigmay]=normfit(y);%计y2均
mu=[mulx muly];%均值向
Sigma=cov(x,y);%协方差矩
[X,Y]=meshgrid(-5:0.1:5,-5:0.1:5); xy=[X(:) Y(:)]; %产生网格数据 p=mvnpdf(xy,mu,Sigma); %求取联合概
Z=reshape(p,size(X));
mesh(X,Y,Z);%画联概率密度三维
title('二维高斯分布的联合概
运行代码之后,可以得
二、画自相关函
因为给定的两个wav文件以采样率22050进行采样,所以在编写程序码之先对两个wav文件的样进行了分析。分的结果得到,high.wav的样正好22050,low.wav的采样率为44100,为22050的倍。也就是说对于high不需要进行重新采样,而对于low也只需要进隔点取,从而降低采样率。然后计算出自相关函
程序代码如下:
clear
[filename,filepath]=uigetfile('.wav','Open wav file');%确定wav文件位置 [y,fs,bits]=wavread([filepath,filename]);%解读wav文件
l=length(y);
if fs==44100 %如果样频率为44100,先抽取偶
a=1:l;
b=a(find(rem(a(:),2)==0));%偶数
x=y(b);
else fs==22050 %果抽样频率为22050,则保持
x=y;
end
z=xcorr(x,'unbiased'); %求自相
subplot(211)
plot(x);
subplot(212)
plot(z);
运行代码之后,可以得
High的波形及其自
Low的波形及其自相关函数曲线 三、画自信号的自相函数曲
这次我用windows自带的录音机录一小段wav文件,并44100的采样进行了分
程序代码如下:
clear
[filename,filepath]=uigetfile('.wav','Open wav file');%定wav文件
[x,fs,bits]=wavread([filepath,filename]);%解读wav文件 z=xcorr(x,'unbiased'); %求自
subplot(211)
plot(x);
subplot(212)
plot(z);
运行代码之后,可以得
经过三个自相关函数的试验,可以看出,自相关函数具偶对称性和收敛,并且自
函数的宽度是原信号
自相关函数与偏自相关函数
自相关函数与偏
上一节介绍了随过程的几种模型。实际中单对时间序列的观察难确定其属于哪一种模型,而自相关数和偏自相关函数是分析随机程识别模型的有力
1、自相关函数定义
在给出自相关函数定之前先介绍自协方差函数概念。由第节知随机过程{x t }中的一个元素x t ,t = 1, 2, … 都是随机变量。对于平稳的随机程,其期望为常数,用μ
E (x t ) =μ,t =1,2,
随机过程的取值将以 μ 为中心上下变动。平稳机过程的方差是一个
2
t =1,2, V a r (x σt ) =x ,
σx 2用来度量随机过程取值对均值μ的离散程
相隔k 期的两个随机变量x t 与x t -k 的协方差滞后k 期的自协差,定义
γ=C o v (x , x ) =E [(x -μ) (x -μ) ]
k
t t -k
t
t -k
, 1,2, 自协方差序列:γk ,k =0
2
V a rx (t ) =σ为随机过程{x t }的自协方差数。当k = 0
2
σ
a r ()x V a r (x ) = 因为对于一个平
t =t -k x
k = 0
=1。
, 1,2, 以滞后期k 为变量的自
o v (x o v (x 即C t -k , x t ) = C t , x t +k ) ,自相关函数是零对称的,所际研究中只给出
函数的正半部
2、自回归过程的自相关函数 (1)平稳AR(1)过程自相关
φx u AR(1) 过程:x t =1t -1+t ,|φ1| <>
已知E (x t ) =0(why? )。用x t -k 乘上式
x t x t -k =φx x u 1t -1t -k +t x t -k
上式两侧同取期望:γk =φ1γk -1
(t t -k ) =0其中Eux (why? )(由
φ12 ut-k-2 +…,而u t 是白噪音与其t - k 期及以前各
两侧同除 γ0 得:
ρ=φρ=φρ==φρ
k
2
1k -11k -2
k
10
k
因为ρo = 1,所以有ρk =φ1(k ≥0)
对于平稳序列有 | φ1| < 1。所以当="" φ1为正时,自相关函数按指数衰至;当="">
1> φ1 > 0 -1<φ1>φ1><>
图 AR(1) 过
同理,对于φ1 =1和φ1 >1情形即非平稳和强平稳过程的自相函数如下
φ1 = 1.1(强非平稳过) φ1 = 1(随机游走
(2)AR(p ) 过
=x +x ++x +u 用x t -k (k > 0) 同乘平稳的 p 阶
x =x x +x x ++x x +x u 的两侧,得:x t -k t 1t -k t -12t -k t -2p t -k t -p t -k t
+++对上式两
用 γ0分别除上式的两侧得Yule-Walker 方
ρk = φ1 ρk -1 + φ2 ρk -2 + … + φp ρk -p , k > 0
φφ
φγ
φ
φφφ
φγφγ
(L ) =1-L -L --L =(1- G L ) 令Φ,其中L 为k 的滞后算子,
i = 1, 2, …, p 是特
-1-2-p k k -1k -2k -p
1-GG -G 0G G +G ,也即G 1i -2i -p i =i =1i +2i +p i 。
φφφ∏
2
12
p
p p
i =1
-1
φφ
k
φφφφ
可证:
ρ=A G +A G ++A G (*)
k
11
k
22
k pp
其中A i , i = 1, … ,p 为定常数。(提
)式代入到Yule-Walker 方程中证明) 由(*)式知道会遇到下几种情
① 当G i 为实数时,(*)式中的AG i i 将随着k 的增加而几何衰减至零,称为指数
k
i i ,G ② 当G i 和G j 表示一对共复数时,设G j =a -b i =a +b G j 的极座标形
R ,则G i ,
GR c o s θ+i s i n θ) i =(
G R (c o s θ-i s i n θ) j =
若AR(p )
R <1。那么随着k>1。那么随着k>
k k GR =(c o s k θ+i s i n) k θ i
k k GR c o s k θ-i s i n) k θ j =(
k
自相关函数(*)式中的相应项G i , G j 将按正弦振形式衰
k
注意:实际中的平稳自回归过程的自关函数常是由指数减和正弦衰减两分混合而
③ 从(*)式可以看出,当特征方程的根值远离单位圆时,k 不必很大,自相关函就会衰减
④ 有一个实根接近1时,自相关函将衰减的很慢,似于线性衰减。当有两个以上根取值接近1时,自相关函同样会衰减的
两个特征根为实根 两个特征为共轭
AR(2) 过程
3、移动平均过程的自相关函数 (1)MA(1) 过程的相关函
u θu 对于MA(1)过程x t =t +1t -1,
γ=E (x x )[=+E (u θu ) (u +θu ) ]
k
t t -k
t
1t -1t -k 1t -k -1
当k = 0时,
γ=E (x x ) =E [(u +θu ) (u +θu ) ]=E (u +2θu u +u )
t t 2
t 1t -1t 1t -1
2
t
1t t -1
2
t -1
1+θσ =(1)
2
当k = 1时,
γ=E (x x ) =+E [(u θu ) (u +θu ) ]
1
t t -1
t
1t -1t -11t -2
22
=E (u u +u +u u +u u ) =θ1σ2 t t -11t -11t t -21t -1t -2
θθθ
t
当 k > 1 时,
γ=E (x x )[=+E (u θu ) (u +θu ) ]=E (u u +θu u +θu u +θu u ) =0
k
t t -k
1t -1t -k 1t -k -1
2
t t -k 1t t -k -11t -1t -k 11t -t -k -1
综合以上三种情形,MA(1)过程自相关函
ρk
θ1 > 0 θ1 <>
图 MA(1)过程
可见MA(1) 过程的自相关函数
(2) MA(q ) 过程的自相函数 MA(q ) 过程自相关函
ρk = 当k > q ,ρk = 0,说明 ρk , k = 0, 1, … 具有截尾
例如,对于MA(2) 过,自相关函
ρ1 ρ2 ρk = 0, k > 2。
4、 ARMA (1, 1) 过程的自相关
ARMA (1, 1) 过程的
的符号取决于 (φ1 -θ1 ) 。 φ1 > 0,指数衰减是平滑的,或正或负。 φ1 <>
对于ARMA (p , q ) 过程,p , q ≥ 2时,自相关函数的表现形式比复杂,可能是指数衰减、正弦减或二者的混合
5、相关图(correlogram ,或估计自相关函数,样本相关函数) 对于一个有限间序列(x 1, x 2, …, x T )用样本
估计总体均值 μ,用样本方差
s = 2
估计总体方差σx 2。
当用样本矩估计随机过程的自相关函数,则称其为相图或估计的自相函数,
r k
k = 0, 1 , 2, …, K , ( K < t="" )="">
r k 是对ρk
C k
是对γk 的估计。
k = 0, 1, 2, …, K ,
C 0 是对γ0的估计。T 是时间序数据的样本容量。实中T 不应太小,最
注意:C k 为有偏估计量。但在样本条件下更有
相关图是对自相关函数的估计。由于MA 过程和ARMA 过程中的MA 分量的自函数具有截尾特性,所以通过相关可以估MA 过程的阶数q 。相关图是识别MA 过程阶数和ARMA 过程中MA 分量阶数的一个重要方法。于年度时间列据,相关图一般取k = 15
r k 的方差近
1。所以在观察相关图时,若r k 的绝对值超
个标准差),就被认为是显地不为零。
第五节 偏
偏自相关函数描述随机过程结构特征另一种方法。用 φkj 表示k 阶自回归过程第j 个回归系数,则k 自回归模型表
x =x +x ++x +u t k 11t -kt 22-k k t -k t
其中φkk 是最后一个回归系数。若把φkk 作是滞后期k 函数,
φφφ
,2, φkk ,k =1
为偏自相关函数。它由下
x φx u t =11t -1+1t
x φx +φx +u t =21t -122t -22t
x =x +x ++x +u t k 11t -kt 22-k k t -k k t
因偏自相关函数中每一个回归系φkk 恰表示x t 与x t -k 在排除了其中间变
φφφ
x t -2,,x t -k +1 影响之后的相关系
xx --x --x =x +u t k 1t -1k 2t -2k k -1t -k +1k k t -k k t
所以偏自相关函数由此得名。 用φkj 表达Yule-Walker
j
k 1j -1k 2j -2
φφφφ
k
ρ=φρφ+ρ++φρ, 得
1k -12k -2
p k -p
k kj -k
ρ=φρφ+ρ++φρ
用矩阵形式表
?ρ1?ρ2
??... ??ρk
或
?1?
?ρ?
?= ?1??... ??
ρ?k -1?
... ρ??φk 1k -1
??φ1... ρk -2??k 2??... ... ... ... ... ??
ρρ... 1??φkk k -2k -3
ρ1
ρ2
ρ1
?
?
? ???
ρ = P φ. 则
φ = P -1ρ,
将k = 1, 2 , … 代入式连续求解,可求得
11 = ρ1,
?φ21??1?φ?=??22??ρ
1
?ρ1?ρ1?1??
-1
?ρ??2?
其中
φ 22 = …
对于AR(1)过程,x t = φ11 x t -1 +u t ,当k = 1时,φ11 ≠ 0;k > 1,φkk =0。所以AR(1)过程的偏自相关函数特征是在k = 1现值(φ11 = ρ1)然
φ11 > 0 φ11 <>
AR(1) 过程
对于AR(2)过程,当k ≤ 2,φkk ≠0;当k >2时,φkk =0。偏自相关函数在滞后2以后有截尾
对于AR(p ) 过程,当k ≤ p 时,φkk ≠0;当k > p时,φkk =0。偏自相关函数在滞后期p 以后有截尾特性,因此可用此特征别AR(p ) 过程的
对于MA(1)过程x t =u t + θ1 ut -1,有 [1/ (1+ θ1 L )]x t =u t , (1- θ1 L + θ12 L 2 - … )x t =u t ,
x t = θ1 x t-1 - θ12 x t-2 +θ13 x t-3 - … +u t
当θ1 > 0时,自回归系数的符号是负交替的;当θ1 <>
因为MA(1) 过程可以转换无限阶的AR 过程,所以MA(1) 过程的偏自相关函数指数衰减
θ1 > 0 θ1 <>
MA(1) 过程的
对于MA(2) 程,若Θ (L ) = 0根是实数,偏自相关数由两个指数衰减形式叠加而成。若Θ (L ) = 0的根是复数,自关函数呈正弦衰减
因为任何一个可的MA(q ) 过程都以转换成一个无限的系数按几何递减的AR 过程,以MA(q ) 过程的偏自函数呈缓慢衰减
ARMA( p , q ) 过程的偏自相关函数也是无限延长的,其表现形式与MA(q ) 程的自关函数相类。根据模型中移动平均部分的数q 以及参数θi 的不同,偏自相关函数指衰减和(或)正弦衰减混
对于时间序列数据,偏自相关函数通常是知的。可以用样本计 φ11, φ22, …
??量 φ11, φ22, …。估计的偏自相关
? φkk , k = 1, 2, …, K ,
称为偏相关图。因AR 过程和ARMA 过程中AR 分量的偏自相关函数有截特性,所以可利用偏相关图估计自回过程的阶数p 。实际中对于偏相关取k = 15就足可
-1?φkk 方差近似为T 。当T 充
-1/2?所以在观察偏相关图时,若φ(2个标准差),就被认是显著地不kk 的
为零。
注:2个标准差 = 2 T -1/2 = 2(1/7)= 0.286。图中虚线表示到中心2个标准差
补充知识:检验过程是否为白声的Q 统
在介绍Q 统计量之前,先介绍序列y t 的估计的相关函数(相图)的定
r k
k = 1, 2, ….
其中r k 表示y t 与y t-k 估计的自相关系,是对自相关系数ρk 的估
∑y t -k )/ (T -k )
∑y t )/ T 。
模型残差序列是否为白噪声的检验用Box-Pierce (1970) 出的Q 统计量完成的。Q 检验的零
H :ρ1 = ρ2 = … = ρK = 0
即序列是一个白噪声过程。其中ρi 表示自关系数。Q 计量定
Q = T
∑r
k =1
K
k
2
(r k 是用残差序列计算的自关系数的估计
2
K -p -q ) 分布,其中T 表示样容量,K 表示自相关数的随着T →∞,Q 渐近服
个数,p 表示模型自回归部分的最大滞后值,q 示移动平均部分最大滞后
K -p -q ) 分布存在差异(相应值Ljung 和Box 认为定义的Q 统计
小),于是提出修正的Q 统计量。
Q = T (T
K -p -q ) 分布。且它其中r k ,K ,p ,q 的定如上式。修正的Q 统
近似性比原Q 统计量的近似性更。(注意:EViews 中给出的Q 统计量就是按修正的Q 统计量定义
用残差序列计算Q 统计量的值。显然若残差列不是白噪声,残差列中必含有其他成份,自相关系数不于零。则Q 值将很大,反之Q 值将很小。判别规
χK -p -q ) 若Q ≤,则接受H 0。 α(2
K -p -q ) ,则拒绝H 0。 若Q >χα(2
其中α 表示验水平;p ,q 分表示时间序列模中自回归和移动平均滞后项的个。 实际检验中,K 取15左右
11