matlab计算二维自相关与功率谱（两者的互推）

wqsoooooooooo · 发表于 2009-1-11 14:29

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

[M,N]=size(A);
R=xcorr2(A)/(M*N)/Sq^2; %%% 其中Sq为A的均方根偏差。

这样对否？

再求功率谱密度怎么求？是不是
S1=fft2(R); S=S1*conj(S1);

这样对不？

[ 本帖最后由 ChaChing 于 2009-4-14 19:47 编辑 ]

ch_j1985 · 发表于 2009-1-11 22:33

在版面顶部使用【搜索】功能搜索一下：自相关函数
http://forum.vibunion.com/forum/ ... 0%B9%D8%BA%AF%CA%FD

wqsoooooooooo · 发表于 2009-2-17 17:11

>> A

A =
   1    2
   3    4

>> P=fft2(A,3,3);
>> P=P.*conj(P);
>> P

P =
  100.0000 28.0000 28.0000
37.0000 13.0000 7.0000
37.0000 7.0000 13.0000

>> P=fftshift(P)

P =
13.0000 37.0000 7.0000
28.0000  100.0000 28.0000
7.0000 37.0000 13.0000

>> R=xcorr2(A);
>> R

R =
   4 11    6
14 30 14
   6 11    4

>> Pr=fft2(R);

>> Pr3=Pr.*conj(Pr)

Pr3 =
  1.0e+004 *
1.0000 0.0784 0.0784
0.1369 0.0169 0.0049
0.1369 0.0049 0.0169

>> Pr3=fftshift(Pr3)

Pr3 =
  1.0e+004 *
0.0169 0.1369 0.0049
0.0784 1.0000 0.0784
0.0049 0.1369 0.0169

>> Pr3=sqrt(Pr3)

Pr3 =
13 37    7
28 100 28
   7 37 13

>> p=fft2(A,3,3)

p =
  10.0000          1.0000 - 5.1962i 1.0000 + 5.1962i
  -0.5000 - 6.0622i  -3.5000 - 0.8660i 2.5000 - 0.8660i
  -0.5000 + 6.0622i 2.5000 + 0.8660i  -3.5000 + 0.8660i

>> p=p.*conj(p)

p =
  100.0000 28.0000 28.0000
37.0000 13.0000 7.0000
37.0000 7.0000 13.0000

>> R=abs(real(ifft2(p)))

R =
30.0000 14.0000 14.0000
11.0000 4.0000 6.0000
11.0000 6.0000 4.0000

>> R=fftshift(R)

R =
4.0000 11.0000 6.0000
14.0000 30.0000 14.0000
6.0000 11.0000 4.0000

[ 本帖最后由 ChaChing 于 2009-4-14 19:50 编辑 ]

wqsoooooooooo · 发表于 2009-2-17 21:14

忘了
P是指功率谱。。
R是自相关

giftdreamer · 发表于 2009-2-17 21:33

:@) 正发愁呢，谢谢了！

ChaChing · 发表于 2009-2-17 21:58

可以合并的尽量合并, 别介意!
针对3F的程序, 两点建议
1.感觉太长不易阅读, 是否考量程序/输出两者分开
2.若能够增加一些说明, 或许就更完善

wqsoooooooooo · 发表于 2009-2-24 19:28

A=[1 2 3;1 2 3;1 2 3];
>> P=fft2(A,5,5);    %%5为A的size的二倍减一。
>> P=P.*conj(P);    %%有A直接求出功率谱P；
>> R=real(ifft2(P));
>> R=fftshift(R)       %%有功率谱再求出自相关R；

R =
3.0000 8.0000 14.0000 8.0000 3.0000
6.0000 16.0000 28.0000 16.0000 6.0000
9.0000 24.0000 42.0000 24.0000 9.0000
6.0000 16.0000 28.0000 16.0000 6.0000
3.0000 8.0000 14.0000 8.0000 3.0000

>> R1=xcorr2(A)             %%有A直接求得的自相关R1

R1 =
   3    8 14    8    3
   6 16 28 16    6
   9 24 42 24    9
   6 16 28 16    6
   3    8 14    8    3

>> P1=fft2(R1);                   %%有自相关R1求功率谱P1
>> P1=P1.*conj(P1);
>> P1=sqrt(P1);
>> P1=fftshift(P1)

P1 =
1.1115 5.3820 13.7508 5.3820 1.1115
7.6180 36.8885 94.2492 36.8885 7.6180
26.1885  126.8115  324.0000  126.8115 26.1885
7.6180 36.8885 94.2492 36.8885 7.6180
1.1115 5.3820 13.7508 5.3820 1.1115

>> P=fftshift(P)

P =
1.1115 5.3820 13.7508 5.3820 1.1115
7.6180 36.8885 94.2492 36.8885 7.6180
26.1885  126.8115  324.0000  126.8115 26.1885
7.6180 36.8885 94.2492 36.8885 7.6180
1.1115 5.3820 13.7508 5.3820 1.1115

现在重新整理了一下！！  呵呵希望有帮助！

[ 本帖最后由 ChaChing 于 2009-4-14 19:34 编辑 ]

单刀赴会 · 发表于 2009-4-12 14:17

哈哈，不错！不过小弟在看代码的时候发现其实楼主可以省去第四行代码： R=real(ifft2(P));
因为A为实信号，所以它的傅里叶变换P是偶函数，conj(P)也是偶函数，第3行P=|P|.^2后P既是偶函数又是实函数，既P是实偶频谱，实偶频谱的傅里叶反变换为实偶信号，故ifft2(P)为实的，所以R=real(ifft2(P))没有起到作用。我用matlab验证了，我的猜想是正确的。不过我还没有想通为何楼主这行代码的最初用意是什么，请指点一下。

wqsoooooooooo · 发表于 2009-4-12 14:43

我是为了实现功率谱来的
发现用xcorr做向相关函数速度很慢！用功率谱再求自相关函数就快多了！
我也不知道其中的道理
呵呵

单刀赴会 · 发表于 2009-4-14 16:56

你的意思是：matlab调用xcorr2函数耗时多，而通过求傅里叶变换再求功率谱然后在傅里叶反变换得到互相关函数更快一些，是吗？
楼主是高手，做个朋友吧？我的qq：674411410.有机会的话，私下请教你。

wqsoooooooooo · 发表于 2009-4-14 17:48

高手不算就光知道那么一点点！

ChaChing · 发表于 2009-4-14 20:00

楼上两位都不错嘛!
可以的话给个建议, 不用私下讨论嘛! LZ讨论时, 不介意我们其他人顺道学习学习!

还有LZ有没发现7F的R及R1输出的差异! why?
是否第一种方法的计算过程较繁杂造成时间/精度的差异!? 个人直觉/猜测!
待高手路过!

[ 本帖最后由 ChaChing 于 2009-4-14 20:55 编辑 ]

ChaChing · 发表于 2009-4-14 20:39

原帖由 wqsoooooooooo 于 2009-4-12 14:43 发表

 登录/注册后可看大图

...用xcorr做向相关函数速度很慢！用功率谱再求自相关函数就快多了！...

学理以前就没学好, 也忘了一乾二净了!
用LZ的程序试了下! 好像刚好相反!? why?

N=3; A=[1:N]; A=repmat(A,N,1); NN=2*N-1;
tic
P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P);
toc
tic
R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1);
toc

复制代码

wqsoooooooooo · 发表于 2009-4-15 09:06

这应该和矩阵的大小有关系吧！
我说的那种方法能用来计算300X300以上的矩阵的功率谱和子相关比较快的！
clear; N=3; A=1:N; A=repmat(A,N,1); NN=2*N-1;
tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.002654 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 0.002915 seconds.
矩阵较小时计算时间相差不大！！
下面把N设为100
clear; N=100; A=1:N; A=repmat(A,N,1); NN=2*N-1;
>> tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.042852 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 0.189567 seconds.
下面把N设为300
clear; N=300; A=1:N; A=repmat(A,N,1); NN=2*N-1;
>> tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.319294 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 12.937269 seconds.
这样我们可以看到速度的差别了啊呵呵
不过应该不同的计算机也可能不同吧！

[ 本帖最后由 ChaChing 于 2009-4-15 21:52 编辑 ]

ChaChing · 发表于 2009-4-15 22:00

的确如此, 昨晚我仅随意试到10而已!
>> clear; N=3; A=1:N; A=repmat(A,N,1); NN=2*N-1;
tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.039205 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 0.010360 seconds.
clear; N=10; A=1:N; A=repmat(A,N,1); NN=2*N-1;
>> tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.013720 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 0.001280 seconds.
clear; N=100; A=1:N; A=repmat(A,N,1); NN=2*N-1;
tic; P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
R=fftshift(R); P=fftshift(P); toc
Elapsed time is 0.053816 seconds.
>> tic; R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
P1=sqrt(P1); P1=fftshift(P1); toc
Elapsed time is 0.127533 seconds.

[ 本帖最后由 ChaChing 于 2009-4-15 22:03 编辑 ]

账号		自动登录	找回密码
密码			我要加入

[编程技巧] matlab计算二维自相关与功率谱（两者的互推）

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复楼主 wqsoooooooooo 的帖子

终于把二维的自相关函数与功率谱搞明白了！（两者之间的互换）

评分

回复板凳 wqsoooooooooo 的帖子

二维自相关函数与功率谱的互转

评分

哈哈，不错

评分

回复 8楼单刀赴会的帖子

回复 9楼 wqsoooooooooo 的帖子

回复 10楼单刀赴会的帖子

评分

回复 13楼 ChaChing 的帖子

评分

[编程技巧] matlab计算二维自相关与功率谱 （两者的互推）

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复 楼主 wqsoooooooooo 的帖子

终于把二维的自相关函数与功率谱搞明白了！（两者之间的互换）

评分

回复 板凳 wqsoooooooooo 的帖子

二维自相关函数与功率谱的互转

评分

哈哈，不错

评分

回复 8楼 单刀赴会 的帖子

回复 9楼 wqsoooooooooo 的帖子

回复 10楼 单刀赴会 的帖子

评分

回复 13楼 ChaChing 的帖子

评分

[编程技巧] matlab计算二维自相关与功率谱（两者的互推）

回复楼主 wqsoooooooooo 的帖子

回复板凳 wqsoooooooooo 的帖子

回复 8楼单刀赴会的帖子

回复 10楼单刀赴会的帖子