高手指导

t12211017

老师给了一个小程序让自己编出来不知道对不对，最后要求数据的均方差，结果应该是与 h=au^4/(4bg^2)的结果应该是一致的，a=0.0081,u是自己设定的风速，不.74，g是重力加速度，我求的h=0.017结果算出来总是不一样啊


clear
N1=128;dx=1;N2=128;a=0.0081;b=0.74;g=9.81;u=5;
y1=randn(N1,N2);
y1=y1-mean2(y1);
y1=y1/sqrt(var(y1(:)));
out=fft2(y1);
outr=real(out);
k1=0;
for n1=1:1:N1
k1=2*pi/N1/dx+k1;
k2=0;
  for n2=1:1:N2
      k2=2*pi/N2/dx+k2;
   wpm(n1,n2)=a/8/(k1^2+k2^2)^2*exp(-b*g*g/(k1^2+k2^2)/u^4);
  end
end
for n1=1:1:N1
    for n2=1:1:N2
    t(n1,n2)=wpm(n1,n2)*out(n1,n2);
    end
end
uuu=ifft2(t);
uu=real(uuu);
smooth(uu);
surf(uu)
shading interp
s1=mean2(uu);
s=uu-s1;
nn=s.*conj(s);
o=sum(nn(:));
y=sqrt(o/size(s,1))