cross spsctrum using cpsd of 2D matrices

Hi there

I'm struggling with cpsd for 2D matrices. In my test, I've generated two 100x100 matrices (sine waves) s1 and p1 with a small phase difference between them. I'm trying to identify the propagation direction using this phase difference, hence using cpsd

[Psp F] = cpsd(s1,p1....).

My guess is that cpsd is the one for the job. However I'm not sure if I'm doing the right thing. When the crests of s1 and p1 are orthogonal to the axis, my script works just fine. However, if they have some angle, instead of two components, lets say in the 1st and 3rd quadrant, I end up with 4 components, in all quadrants. Not good...

I suppose there's sth wrong with the way I'm using cpsd for matrices. Anyone out there could give me a hand? Many thanks.

This is my code so far

[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); % 

[XF YF] = meshgrid(F, F);      

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

figure, surface(XF,YF,abs(Psp))

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