computing the cross spectrum of two-dimensional matrice

Hi there

I got stuck with the matlab syntax trying to compute the cross spectrum of two-dimensional matrices.

I suppose that what I need is something like

[Pxy,F]=cpsd(x,y,[],[],[],dt);

where x and y are 100x100 matrices. Obviously the above syntax doesn't work, but that works for my vectors.

I've browsed the web and couldn’t find anything really.

Is there anyone out there kind enough to show me the correct syntax? Or is there any other command that does the job instead of cpsd?

I'm working with MATLAB Version 7.11.0.584 (R2010b)

Tips are highly welcomed, regards

nvc

Answer 1

[Walter Roberson]

That syntax is documented for R2015a. http://www.mathworks.com/help/signal/ref/cpsd.html

<<The input signals may be either vectors or two-dimensional matrices. If both are vectors, they must have the same length. If both are matrices, they must have the same size, and cpsd operates columnwise: Pxy(:,n) = cpsd(x(:,n),y(:,n))>>

It appears that support was added in R2014b. You can use a loop for earlier versions.

[n vc]

Dear Roberson

thank you for your response.

i have tried sth like

[l c]=size(s1);

  for c=1:c

[Psp1(:,c),F]= cpsd(s1(:,c),p1(:,c),[],[],[],dt);

  end

  for l=1:l

[Psp2(l,:),F]= cpsd(s1(l,:),p1(l,:),[],[],[],dt);

  end

Psp=Psp1*Psp2;

for the 1st quadrant it seems to work, but I need to compute the negative frequencies as well. My main goal is to detect the phase difference between signals to eliminate the directional ambiguity (moving from or towards?).

Can you give me some advice?

regards

[Walter Roberson]

for c=1:c and for l=1:l might work in practice but they are going to confuse the heck out of people reading them. You should avoid using a loop control variable named the same thing as one of the bounds of the loop. For example,

[numrow numcol] = size(s1);
for c = 1:numcol
  [Psp1(:,c),F]= cpsd(s1(:,c),p1(:,c),[],[],[],dt);
end
for l = 1:numrow
  [Psp2(l,:),F]= cpsd(s1(l,:),p1(l,:),[],[],[],dt);
end

For the other quadrants: does the 'twosided' parameter help with that?

[n vc]

Dear Roberson

yes, twosided seems to be the way, this is how it looks now:

[numrow numcol] = size(s1);

  for c = 1:numcol

[Psp1(:,c),F]= cpsd(s1(:,c),p1(:,c),[],[],[],dt,'twosided');

  end

  for l = 1:numrow

[Psp2(l,:),F]= cpsd(s1(l,:),p1(l,:),[],[],[],dt,'twosided');

  end

Psp=Psp1*Psp2;

F1=F-1/(2*dt); 

[XF1 YF1] = meshgrid(F1, F1);            

figure, surface(XF1,YF1,fftshift(abs(Psp)))

However, there's still sth wrong. My component should be only in the 1st and 3rd quadrants, but it appears as well on the 2nd and 4th, not clear why. Any more ideas on that?

Thanks again for your time, it's been a valuable contribution.