互功率谱与互相关函数的一个编程问题

yangcui

以下是我作的，遇到困难了，请大家帮忙解决一下。求两个信号序列的互相关函数，4.txt和6.txt中第二列是实测信号序列，第一列为采样时间。我按照一本书所说的，先求出两者的互功率谱，然后进行傅立叶逆变换，取傅立叶变换的实部为互相关函数，结果如下图所示。我感觉和预想的结果不一样。下面是用的程序，是不是建立负频率段的互功率谱错了，还是其他的地方出错了，请高手帮忙。万分感激！急用！！！
clc
clear
data1=load('4.txt');
data2=load('6.txt');
x=data1(:,2)';
y=data2(:,2)';
np=length(x);
fs=1000;
nfft=1024;
t=0:1/fs:(np-1)/fs;
Pxy=CSD(x,y,nfft,fs,[]);  %计算互功率谱
Pxy(nfft/2+1)=real(Pxy(nfft/2));  %建立负频率段的互功率谱
Pxy(nfft/2+2:nfft)=conj(Pxy(nfft/2:-1:2));
g=ifft(Pxy);    %傅立叶逆变换
g=ifft(Pxy);    %取傅立叶变换的实部为互相关函数
r=real(g(1:np));%绘制互相关函数时程曲线图
plot(t,r);
grid;

[ 本帖最后由 zhangnan3509 于 2007-6-5 17:19 编辑 ]

yangcui

上面的程序贴错了，应该是下面的这个。
clc
clear
data1=load('4.txt');
data2=load('6.txt');
x=data1(:,2)';
y=data2(:,2)';
np=length(x);
fs=1000;
nfft=1024;
t=0:1/fs:(np-1)/fs;

Pxy=CSD(x,y,nfft,fs,[]);  %计算互功率谱
Pxy(nfft/2+1)=real(Pxy(nfft/2));  %建立负频率段的互功率谱
Pxy(nfft/2+2:nfft)=conj(Pxy(nfft/2:-1:2));
g=ifft(Pxy);    %傅立叶逆变换
g=real(ifft(Pxy));    %取傅立叶变换的实部为互相关函数
r=g(1:np);%绘制互相关函数时程曲线图
plot(t,r);
grid;

songzy41

x和y的互相关函数应有2001个数据，很明显没法从CSD得到。

yangcui

songzy41，您好！能具体指点一下吗？
在什么条件下能从CSD得到？
我如果按下面的步骤作，能行吗？
1、将x，y长度通过增加零值延长到2001；
2、利用FFT计算x和y各自2001点的DFT：X（k）和Y（k）；
3、计算R（k）＝X（k）Y‘（k），其中Y‘（k）为Y（k）的共轭；
4、利用IFFT计算R（k）得到x，y的互相关函数。

yangcui

songzy41，您好！能具体指点一下吗？
在什么条件下能从CSD得到？
我如果按下面的步骤作，能行吗？
1、将x，y长度通过增加零值延长到2001；
2、利用FFT计算x和y各自2001点的DFT：X（k）和Y（k）；
3、计算R（k）＝X（k）Y‘（k），其中Y‘（k）为Y（k）的共轭；
4、利用IFFT计算R（k）得到x，y的互相关函数。

songzy41

原帖由 yangcui 于 2007-1-29 22:51 发表
1、将x，y长度通过增加零值延长到2001；
2、利用FFT计算x和y各自2001点的DFT：X（k）和Y（k）；
3、计算R（k）＝X（k）Y‘（k），其中Y‘（k）为Y（k）的共轭；
4、利用IFFT计算R（k）得到x，y的互相关函数。

用这种方法可以求得互相关函数，但补零时要注意，x应在后面补零，y应在前面补零。程序如下：
clc; clear
data1=load('4.txt');
data2=load('6.txt');
x=data1(:,2)';
y=data2(:,2)';
np=length(x);
x=[x zeros(1,1000)];
y=[zeros(1,1000) y];
fs=1000;
X=fft(x);
Y=fft(y);
R=X.*conj(Y);
rxy=real(ifft(R));
t=(-np+1:np-1)/fs;
plot(t,rxy); grid;
得图有：

吴向南

x和y的互相关函数应有2001个数据，很明显没法从CSD得到。
请问x和y的互相关函数为什麽是2001个数据呢？而不是1001个数据呢？？本人不是学信号的，但是现在要用到这个程序，学习下，问题比较初级，请帮忙回答一下，不甚感激。


补充内容 (2012-10-13 21:08):
还用csd函数，把楼主的程序改成如下：对不对呢？？？请赐教~~~~~不甚感激
clear
data1=load('4.txt');
data2=load('6.txt');
x=data1(:,2)';
y=data2(:,2)';
np=length(x);
fs=1000;
nfft=2^nextpow2(2*np);
t=0:1/fs:(np-1)/fs;
Pxy=CSD(x,y,nfft,fs,[]);  %计算互功率谱
Pxy(nfft/2+1)=real(Pxy(nfft/2));  %建立负频率段的互功率谱
Pxy(nfft/2+2:nfft)=conj(Pxy(nfft/2:-1:2));
g=ifft(Pxy);    %傅立叶逆变换
g=real(g(1:np));    %取傅立叶变换的实部为互相关函数
r=g(1:np);
plot(t,r);
grid;

anpengqian

求指点。。。这两句具体怎么理解？
Pxy(nfft/2+1)=real(Pxy(nfft/2));  %建立负频率段的互功率谱
Pxy(nfft/2+2:nfft)=conj(Pxy(nfft/2:-1:2));