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 编辑 ]
回复 楼主 wqsoooooooooo 的帖子
在版面顶部使用【搜索】功能搜索一下：自相关函数http://forum.vibunion.com/forum/viewthread.php?tid=36108&highlight=%D7%D4%CF%E0%B9%D8%BA%AF%CA%FD
终于把二维的自相关函数与功率谱搞明白了！（两者之间的互换）
>> AA =
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 编辑 ] 忘了
P是指功率谱。。
R是自相关 :@) 正发愁呢，谢谢了！
回复 板凳 wqsoooooooooo 的帖子
可以合并的尽量合并, 别介意!针对3F的程序, 两点建议
1.感觉太长不易阅读, 是否考量程序/输出两者分开
2.若能够增加一些说明, 或许就更完善
二维自相关函数与功率谱的互转
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 编辑 ]
哈哈，不错
哈哈，不错！不过小弟在看代码的时候发现其实楼主可以省去第四行代码： R=real(ifft2(P));因为A为实信号，所以它的傅里叶变换P是偶函数，conj(P)也是偶函数，第3行P=|P|.^2后P既是偶函数又是实函数，既P是实偶频谱，实偶频谱的傅里叶反变换为实偶信号，故ifft2(P)为实的，所以R=real(ifft2(P))没有起到作用。我用matlab验证了，我的猜想是正确的。不过我还没有想通为何楼主这行代码的最初用意是什么，请指点一下。
回复 8楼 单刀赴会 的帖子
我是为了实现功率谱来的发现 用xcorr做向相关函数速度很慢！ 用功率谱再求自相关函数就快多了！
我也不知道其中的道理
呵呵
回复 9楼 wqsoooooooooo 的帖子
你的意思是：matlab调用xcorr2函数耗时多，而通过求傅里叶变换再求功率谱然后在傅里叶反变换得到互相关函数更快一些，是吗？楼主是高手，做个朋友吧？我的qq：674411410.有机会的话，私下请教你。
回复 10楼 单刀赴会 的帖子
高手不算就光知道那么一点点！ 楼上两位都不错嘛!可以的话给个建议, 不用私下讨论嘛! LZ讨论时, 不介意我们其他人顺道学习学习!
还有LZ有没发现7F的R及R1输出的差异! why?
是否第一种方法的计算过程较繁杂造成时间/精度的差异!? 个人直觉/猜测!
待高手路过!
[ 本帖最后由 ChaChing 于 2009-4-14 20:55 编辑 ] 原帖由 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
回复 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 编辑 ] 的确如此, 昨晚我仅随意试到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