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

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

=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

回复 楼主 wqsoooooooooo 的帖子

在版面顶部使用【搜索】功能搜索一下:自相关函数
http://forum.vibunion.com/forum/viewthread.php?tid=36108&highlight=%D7%D4%CF%E0%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.0000100.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

回复 板凳 wqsoooooooooo 的帖子

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

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

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

A=;
>> 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.1885126.8115324.0000126.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.1885126.8115324.0000126.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

回复 8楼 单刀赴会 的帖子

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

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

回复 9楼 wqsoooooooooo 的帖子

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

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

回复 10楼 单刀赴会 的帖子

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

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 发表 http://www.chinavib.com/forum/images/common/back.gif
...用xcorr做向相关函数速度很慢! 用功率谱再求自相关函数就快多了!...
学理以前就没学好, 也忘了一乾二净了!
用LZ的程序试了下! 好像刚好相反!? why?
N=3; A=; 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

回复 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 编辑 ]

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 编辑 ]
页: [1] 2
查看完整版本: matlab计算二维自相关与功率谱 (两者的互推)