fp=fopen（‘delg_new.grd‘‘r‘）;
fscanf（fp‘%s‘1）;
nx=fscanf（fp‘%d‘1）;
ny=fscanf（fp‘%d‘1）;
xmin=fscanf（fp‘%g‘1）;
xmax=fscanf（fp‘%g‘1）;
ymin=fscanf（fp‘%g‘1）;
ymax=fscanf（fp‘%g‘1）;
fscanf（fp‘%g‘2）;
for?i=1:ny
????for?j=1:nx
????????data（ij）=fscanf（fp‘%g‘1）;
????end
end
fclose（fp）;

dx=（xmax-xmin）/（nx-1）;
dy=（ymax-ymin）/（ny-1）;
x=xmin:dx:（xmin+（nx-1）*dx）;
y=ymin:dy:（ymin+（ny-1）*dy）;

[X?Y]=meshgrid（xy）;
%figure（1）contour（XYdata）xlabel（‘X（m）‘）ylabel（‘Y（m）‘）title（‘原始重力數據‘）;
figure（1）surf（Xdata）view（090）;colorbar;shading?interp;

[CS]=wavedec2（data2‘haar‘）;
Yh2=detcoef2（‘h‘CS2）;
%?Yv2=detcoef2（‘v‘CS2）;
%?Yd2=detcoef2（‘d‘CS2）;
figure（2）surf（X/2Yh2）view（090）;colorbar;shading?interp;
%?figure（3）contour（XYYv2）xlabel（‘X（m）‘）