function z=q(x,y)
z=0;m=200;W=2;n0=0.1;n1=0.011;n2=2.83;Gq0=256*1e-6;
for i=1:m % m为频率范围【n1,n2】划分的m个小区间
n=n1+(i-0.5)*(n2-n1)/m; % n为每个小区间内中心频率
Gq=Gq0*(n/n0)^(-W);
z=z+sqrt(2*Gq*(n2-n1)/m)*sin(2*pi*n*(x+y)+2*pi*rand()); % Gq0为Gq(n0) Gq=Gq0*(n/n0)^(-W) W=2,n0=0.1
end
end