声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6460|回复: 15

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

[复制链接]
发表于 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 编辑 ]
回复
分享到:

使用道具 举报

发表于 2009-1-11 22:33 | 显示全部楼层

回复 楼主 wqsoooooooooo 的帖子

在版面顶部使用【搜索】功能搜索一下:自相关函数
http://forum.vibunion.com/forum/ ... 0%B9%D8%BA%AF%CA%FD
 楼主| 发表于 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 编辑 ]

评分

1

查看全部评分

 楼主| 发表于 2009-2-17 21:14 | 显示全部楼层
忘了
P是指功率谱。。
R是自相关
发表于 2009-2-17 21:33 | 显示全部楼层
:@) 正发愁呢,谢谢了!
发表于 2009-2-17 21:58 | 显示全部楼层

回复 板凳 wqsoooooooooo 的帖子

可以合并的尽量合并, 别介意!
针对3F的程序, 两点建议
1.感觉太长不易阅读, 是否考量程序/输出两者分开
2.若能够增加一些说明, 或许就更完善
 楼主| 发表于 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 编辑 ]

评分

2

查看全部评分

发表于 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验证了,我的猜想是正确的。不过我还没有想通为何楼主这行代码的最初用意是什么,请指点一下。

评分

2

查看全部评分

 楼主| 发表于 2009-4-12 14:43 | 显示全部楼层

回复 8楼 单刀赴会 的帖子

我是为了实现功率谱来的
发现 用xcorr做向相关函数速度很慢! 用功率谱再求自相关函数就快多了!
我也不知道其中的道理  
呵呵
发表于 2009-4-14 16:56 | 显示全部楼层

回复 9楼 wqsoooooooooo 的帖子

你的意思是:matlab调用xcorr2函数耗时多,而通过求傅里叶变换再求功率谱然后在傅里叶反变换得到互相关函数更快一些,是吗?
楼主是高手,做个朋友吧?我的qq:674411410.有机会的话,私下请教你。
 楼主| 发表于 2009-4-14 17:48 | 显示全部楼层

回复 10楼 单刀赴会 的帖子

高手不算  就光知道那么一点点!
发表于 2009-4-14 20:00 | 显示全部楼层
楼上两位都不错嘛!
可以的话给个建议, 不用私下讨论嘛! LZ讨论时, 不介意我们其他人顺道学习学习!

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

[ 本帖最后由 ChaChing 于 2009-4-14 20:55 编辑 ]
发表于 2009-4-14 20:39 | 显示全部楼层

学理以前就没学好, 也忘了一乾二净了!
用LZ的程序试了下! 好像刚好相反!? why?
  1. N=3; A=[1:N]; A=repmat(A,N,1); NN=2*N-1;
  2. tic
  3. P=fft2(A,NN,NN); P=P.*conj(P); R=real(ifft2(P));
  4. R=fftshift(R); P=fftshift(P);
  5. toc

  6. tic
  7. R1=xcorr2(A); P1=fft2(R1); P1=P1.*conj(P1);
  8. P1=sqrt(P1); P1=fftshift(P1);
  9. toc
复制代码

评分

1

查看全部评分

 楼主| 发表于 2009-4-15 09:06 | 显示全部楼层

回复 13楼 ChaChing 的帖子

这应该和矩阵的大小有关系吧!
我说的那种方法能用来计算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 编辑 ]

评分

1

查看全部评分

发表于 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 编辑 ]
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 11:56 , Processed in 0.099064 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表