[已解决] Matlab求两个序列的互相关函数

[复制链接]
王瀚 发表于 2009-5-27 09:42:11
见代码!
计算两序列的互相关函数
%x(n)=[1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1]
%y(n)=x(n-3)+e(n),e(n)为随机噪声
x=[1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
y=[0 0 0 1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
k=length(y);
e=randn(1,k);
y=y+e;
xk=fft(x,2*k);
yk=fft(y,2*k);
rm=real(ifft(conj(xk).*yk));%作用?
rm=[rm(k+2:2*k) rm(1:k)];%作用??
m=(-k+1):(k-1);%作用??
stem(m,rm)
xlabel('m');
ylabel('幅度');
title('x与y的互相关函数');
grid on
hold on
%end
图我贴出来了,可是怎么分析哪?!比较菜,大侠们见笑了!谢谢回复!

[ 本帖最后由 mooni 于 2009-5-27 10:05 编辑 ]

图.fig

5.07 KB, 下载次数: 94778

最佳答案


songzy41 发表于 2009-5-27 11:23:49
原帖由 王瀚 于 2009-5-27 09:42 发表
见代码!
rm=real(ifft(conj(xk).*yk));%作用?
rm=[rm(k+2:2*k) rm(1:k)];%作用??
m=(-k+1):(k-1);%作用??

这是用FFT的方法求互相关函数。其中的机理不是很容易说明白,建议LZ看一下书,在“数字时间序列分析一书”中有一节“利用快速傅里叶变换的相关函数”,详细说明了这种方法。
首先要知道为什么在FFT变换时就补了0,如果不补0的话也可以求相关函数,但求出的是循环相关函数。
rm=real(ifft(conj(xk).*yk)); 求了互功率谱,再反变换就是互相关系数
rm=[rm(k+2:2*k) rm(1:k)]; 由于补0的安排,要把序列以k+1为中心点,左右换置,相当于执行ifftshift函数。但是我发现这语句中把rm(k+1)丢失了,似乎有点问题
m=(-k+1):(k-1);  互相关函数对应的延迟量

我把程序改一下,得下图:
x=[1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
y=[0 0 0 1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
k=length(y);
e=randn(1,k);
y=y+e;
xk=fft(x,2*k);
yk=fft(y,2*k);
rm=real(ifft(conj(xk).*yk));%作用?
%rm=[rm(k+2:2*k) rm(1:k)];%作用??
rm=ifftshift(rm);
m=-k:k-1;%作用??
stem(m,rm)
xlabel('m');
ylabel('幅度');
title('x与y的互相关函数');
回复此楼

60 条回复


songzy41 发表于 2009-5-27 11:23:49
原帖由 王瀚 于 2009-5-27 09:42 发表
见代码!
rm=real(ifft(conj(xk).*yk));%作用?
rm=[rm(k+2:2*k) rm(1:k)];%作用??
m=(-k+1):(k-1);%作用??

这是用FFT的方法求互相关函数。其中的机理不是很容易说明白,建议LZ看一下书,在“数字时间序列分析一书”中有一节“利用快速傅里叶变换的相关函数”,详细说明了这种方法。
首先要知道为什么在FFT变换时就补了0,如果不补0的话也可以求相关函数,但求出的是循环相关函数。
rm=real(ifft(conj(xk).*yk)); 求了互功率谱,再反变换就是互相关系数
rm=[rm(k+2:2*k) rm(1:k)]; 由于补0的安排,要把序列以k+1为中心点,左右换置,相当于执行ifftshift函数。但是我发现这语句中把rm(k+1)丢失了,似乎有点问题
m=(-k+1):(k-1);  互相关函数对应的延迟量

我把程序改一下,得下图:
x=[1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
y=[0 0 0 1 3 2 6 2 1 -2 0 1 5 3 2 -3 0 1 2 0 3 1];
k=length(y);
e=randn(1,k);
y=y+e;
xk=fft(x,2*k);
yk=fft(y,2*k);
rm=real(ifft(conj(xk).*yk));%作用?
%rm=[rm(k+2:2*k) rm(1:k)];%作用??
rm=ifftshift(rm);
m=-k:k-1;%作用??
stem(m,rm)
xlabel('m');
ylabel('幅度');
title('x与y的互相关函数');
we13a.jpg
回复此楼

upcliujie 发表于 2009-8-28 10:42:26

谢谢

看了楼上的回帖,感觉很有帮助

emil0203 发表于 2009-9-21 13:39:37
对我也挺有帮助的,谢谢!

suiying1983654 发表于 2009-9-22 19:21:08
2#真是个专业人士呀,小弟有空得向你多多请教。

zzpanda5 发表于 2009-9-25 09:29:09
非常感谢songzy41,很有帮助!

chunjiaohu 发表于 2010-4-15 09:28:17
受益耶~~~~~~~~~~~~~~~~~~

bianqiyong 发表于 2010-4-17 14:02:10
谢谢二楼详解

chengp205 发表于 2010-4-28 16:26:39
受益了,非常感谢!

whj-winner 发表于 2010-4-28 16:48:32
:$ 佩服:funk: 实在佩服!学习啦!

小白兔1 发表于 2010-5-1 14:11:54
受益匪浅,谢谢

swingdu 发表于 2010-5-25 15:49:22
收益不浅,严重感谢!

gilaroney 发表于 2012-2-24 23:31:43
:handshake
谢谢~~~

eagle1208 发表于 2012-2-26 22:21:46
handshake

jackxiangjiang 发表于 2012-2-27 11:20:07
很好~,ifftshift  好用

amandavon 发表于 2012-4-8 10:12:27
本帖最后由 amandavon 于 2012-4-8 10:20 编辑
songzy41 发表于 2009-5-27 11:23
这是用FFT的方法求互相关函数。其中的机理不是很容易说明白,建议LZ看一下书,在“数字时间序列分析一书” ...


请问一下,如说是两个等长的时间序列要怎么画出这张图呢?
例如,附件中的future和stock两个时间序列是等长的,求其互现关系数?最近一直很疑惑,请指教!

data.xls

112 KB, 下载次数: 3364


songzy41 发表于 2012-4-8 17:37:54
LZ现计算出的的相关系数如同一个三角形,这完全是直流分量太大有关。改为
u=xlsread('data.xls');
u(:,1)=u(:,1)-mean(u(:,1));
u(:,2)=u(:,2)-mean(u(:,2));
[c,LAGS]=xcorr(u(:,1),u(:,2),'coeff');
plot(LAGS,c);
得的图就对了。

chybeyond 发表于 2012-9-8 21:29:54
二楼牛人,学习了

baiwei19890914 发表于 2013-3-4 20:55:33
谢谢啦,很有帮助的.

Erin13 发表于 2014-4-7 19:18:11
songzy41 发表于 2009-5-27 11:23
这是用FFT的方法求互相关函数。其中的机理不是很容易说明白,建议LZ看一下书,在“数字时间序列分析一书” ...

我想问一下rm是什么含义,e=randn(1,k);y=y+e;这两句是什么意思

songzy41 发表于 2014-4-8 09:24:11
Erin13 发表于 2014-4-7 19:18
我想问一下rm是什么含义,e=randn(1,k);y=y+e;这两句是什么意思

rm是互相关函数,e=randn(1,k);产生随机噪声,y=y+e在信号上叠加随机噪声。

Erin13 发表于 2014-4-9 17:29:09
songzy41 发表于 2014-4-8 09:24
rm是互相关函数,e=randn(1,k);产生随机噪声,y=y+e在信号上叠加随机噪声。

clear all
N=256;
n=2*pi/N;
i=n:n:2*pi;
x=5*sin(i);
psi=0.25*pi;
y=10*sin(i+psi);
k=length(y);
e=randn(1,k);%生成1*k的随机项矩阵
y=y+e;%添加白噪声信号
xk=fft(x,2*k);
yk=fft(y,2*k);%2*k是傅里叶变换的点数
rm=real(ifft(conj(xk).*yk));
rm=ifftshift(rm);%求反变换
m=-k:k-1;
stem(m,rm)%绘制针状图
xlabel('m');
ylabel('幅度');
title('x与y的互相关');
这个程序理论上最大值应该出现在32的位置,为什么不行呢?

songzy41 发表于 2014-4-10 09:04:48
Erin13 发表于 2014-4-9 17:29
clear all
N=256;
n=2*pi/N;

现在得到的最大值在负延迟数值上,把
rm=real(ifft(conj(xk).*yk));
改为
rm=real(ifft(xk.*conj(yk)));
对应的延迟就在正延迟数值上了。

zhxw1991 发表于 2014-7-2 23:15:28
有所收获

smile青茶 发表于 2014-9-2 15:12:49
songzy41 发表于 2014-4-10 09:04
现在得到的最大值在负延迟数值上,把
rm=real(ifft(conj(xk).*yk));
改为

你好,作为菜鸟请教几个问题哈,首先y(n)=x(n-3)+e(n)为什么序列y(n)就是在x(n)的基础上左面添了3个零呢?如果给出n的范围呢?例如:x(n)=【3,9,-2】,-3≤n≤3,y(n)=x(n-2)+w(n),w(n)是均值为0,方差为1的高斯随机序列,y(n)又该怎么表示呢?还有求相关系数的话直接用a=xcorr(x,y)求可不可以呢?谢谢
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表