请教下面的程序，反算结果为何与原结果不同？

heaventian

如果 A=[tempa1,tempa2;tempb1,tempb2];b=[H;W],x=[HI,WI];then b=A*x
再又x=inv（A）*b；
但是我的程序反算的结果却完全错误。求解时也没有出现行列式奇异的警告。
请教大家问题出在了哪里？
非常感谢！

clear
load h0
rhoc=2.67e3;
rhom=3.1e3;
drho=rhom-rhoc;
g=9.8;
E=7e10;
niu=0.25;
Te=50e3;
dx=50e3;     
dy=50e3;

% % % % % % % % % % %  
mean_h0=mean(mean(h0));
h0=h0-mean_h0;
HI=fftshift(fft2(h0));
size_HI=size(HI);
tempN=size_HI(1);
tempM=size_HI(2);
kx=(0:tempM-1)/(tempM-1)/dx-1/2/dx;% kx=kx*2*pi
kx=kx*2*pi;
ky=(0:tempN-1)/(tempN-1)/dy-1/2/dy;% kx=kx*2*pi
ky=ky*2*pi;

temp1=sqrt(kx.^2+ky.^2);
k1=temp1;
[kx,ky]=meshgrid(kx,ky);
k=sqrt(kx.^2+ky.^2);
D=E*Te^3/12/(1-niu^2)*1e-16;
k_temp1=k*1e4;
ksi=1+D*k_temp1.^4/drho/g;
phi=1+D*k_temp1.^4/rhoc/g;

% % % % % % % % % % % %
f=0;
WI=f*rhoc*HI/drho*exp(i*randn(41)*pi);
% % % A=[tempa1,tempa2;tempb1,tempb2];b=[H;W],x=[HI,WI];then b=A*x
tempa1=drho*ksi./(rhoc+drho*ksi);
tempa2=drho./(drho+rhoc*phi);
H=HI.*tempa1+WI.*tempa2;
tempb1=rhoc./(rhoc+drho*ksi);
tempb2=rhoc*phi./(drho+rhoc*phi);  
W=HI.*tempb1+WI.*tempb2;
h=ifft2(ifftshift(H));
w=ifft2(ifftshift(W));

     
temp_delta=tempa1.*tempb2-tempa2.*tempb1;
HI=(H.*tempb2-W.*tempa2)./temp_delta;
WI=(-H.*tempb1+W.*tempa1)./temp_delta;
HT=tempa1.*HI;
HB=tempa2.*WI;
WT=tempb1.*HI;
WB=tempb2.*WI;
ht=ifft2(ifftshift(HT));
hb=ifft2(ifftshift(HB));
wt=ifft2(ifftshift(WT));
wb=ifft2(ifftshift(WB));

如果结果正确，HI应该等于原始的HI，WI应该为零。尤其是hb,wb期待为0，ht，wt期待为有限值。单结果却并非如此。程序中为。mat文件，但是附件中为ascii文件。

heaventian

我将kx,ky改了一下，成为
kx=(0:tempM-1)/tempM/dx-1/2/dx;% kx=kx*2*pi
ky=(0:tempN-1)/tempN/dy-1/2/dy;% kx=kx*2*pi
运行出来的ht等就变为了实数，而非（NAN+i*某个数）
但是，对于fftshift,理论上
kx,ky应该如程序中的样子吧？
if mod(tempN,2)==0
    kx=(0:tempM-1)/tempM/dx-tempM/2/tempM/dx;% kx=kx*2*pi
else
    kx=(0:tempM-1)/tempM/dx-(tempM-1)/2/tempM/dx;% kx=kx*2*pi
end
kx=kx*2*pi;
if mod(tempM,2)==0
    ky=(0:tempN-1)/tempN/dy-tempN/2/tempN/dy;% kx=kx*2*pi
else
    ky=(0:tempN-1)/tempN/dy-(tempN-1)/2/tempN/dy;% kx=kx*2*pi
end
ky=ky*2*pi;
然而此时，结果又是错误的。
运行出来的ht等就变为了（NAN+i*某个数）
请教怎么修改？

clear
load h0
rhoc=2.67e3;
rhom=3.1e3;
drho=rhom-rhoc;
g=9.8;
E=7e10;
niu=0.25;
Te=50e3;
dx=50e3;    %the grid space is 4 minute in longitude(or x direction)
dy=50e3;

% % % % % % % % % % %  
mean_h0=mean(mean(h0));
h0=h0-mean_h0;
HI=fftshift(fft2(h0));
size_HI=size(HI);
tempN=size_HI(1);
tempM=size_HI(2);
if mod(tempN,2)==0
    kx=(0:tempM-1)/tempM/dx-tempM/2/tempM/dx;% kx=kx*2*pi
else
    kx=(0:tempM-1)/tempM/dx-(tempM-1)/2/tempM/dx;% kx=kx*2*pi
end
kx=kx*2*pi;
if mod(tempM,2)==0
    ky=(0:tempN-1)/tempN/dy-tempN/2/tempN/dy;% kx=kx*2*pi
else
    ky=(0:tempN-1)/tempN/dy-(tempN-1)/2/tempN/dy;% kx=kx*2*pi
end
ky=ky*2*pi;

temp1=sqrt(kx.^2+ky.^2);
k1=temp1;
[kx,ky]=meshgrid(kx,ky);
k=sqrt(kx.^2+ky.^2);
D=E*Te^3/12/(1-niu^2)*1e-16;
k_temp1=k*1e4;
ksi=1+D*k_temp1.^4/drho/g;
phi=1+D*k_temp1.^4/rhoc/g;

% % % % % % % % % % % %
f=0;
WI=f*rhoc*HI/drho*exp(i*randn(41)*pi);
% % % A=[tempa1,tempa2;tempb1,tempb2];b=[H;W],x=[HI,WI];then b=A*x
tempa1=drho*ksi./(rhoc+drho*ksi);
tempa2=drho./(drho+rhoc*phi);
H=HI.*tempa1+WI.*tempa2;
tempb1=rhoc./(rhoc+drho*ksi);
tempb2=rhoc*phi./(drho+rhoc*phi);  
W=HI.*tempb1+WI.*tempb2;
% h=real(ifft2(ifftshift(H)))+mean_h0;
% w=real(ifft2(ifftshift(W)))+40e3;
h=ifft2(ifftshift(H));
w=ifft2(ifftshift(W));
save h_w h w
     

     
temp_delta=tempa1.*tempb2-tempa2.*tempb1;
HI=(H.*tempb2-W.*tempa2)./temp_delta;
WI=(-H.*tempb1+W.*tempa1)./temp_delta;
HT=tempa1.*HI;
HB=tempa2.*WI;
WT=tempb1.*HI;
WB=tempb2.*WI;
ht=ifft2(ifftshift(HT));
hb=ifft2(ifftshift(HB));
wt=ifft2(ifftshift(WT));
wb=ifft2(ifftshift(WB));