clc;clear;close all;
g=dlmread('a2.dat'); %导入数据
%基本参数
g=g(:,3);%将导入的数据赋值给g
xn=59; %点数
yn=71; %线数
dx=5; %点距
dy=5; %线距
exn=10; %扩边点数
eyn=10; %扩边点数
bk=mean(g(:)); %背景场
nxn=xn+2*exn; %扩边后点数
nyn=yn+2*eyn; %扩边后线数
nxn2=floor(nxn/2)+1; %半数
nyn2=floor(nyn/2)+1; %半数
dfx=1/(dx*nxn); %x方向基频
dfy=1/(dy*nyn); %y方向基频
%将数据重排为矩阵
g=reshape(g,xn,yn)'-bk;
%扩边
gex=zeros(yn+2*eyn,xn+2*exn);
gex(eyn+1:end-eyn,exn+1:end-exn)=g;
gex=[cos([eyn:-1:1]'/eyn*pi/2)*cos([exn:-1:1]/exn*pi/2)*g(1,1) cos([eyn:-1:1]'/eyn*pi/2)*g(1,:) cos([eyn:-1:1]'/eyn*pi/2)*cos([1:exn]/exn*pi/2)*g(1,end)
g(:,1)*cos([exn:-1:1]/exn*pi/2) g g(:,end)*cos([1:exn]/exn*pi/2)
cos([1:eyn]'/eyn*pi/2)*cos([exn:-1:1]/exn*pi/2)*g(end,1) cos([1:eyn]'/eyn*pi/2)*g(end,:) cos([1:eyn]'/eyn*pi/2)*cos([1:exn]/exn*pi/2)*g(end,end)];
gex=gex+bk;
%计算频谱
Sgex=fft2(gex);
%计算延拓因子
[u,v]=meshgrid([0:(nxn2-1) -floor((nxn-1)/2):-1]*dfx,(0:(nyn2-1))*dfy);
f=exp(2*pi*(u.^2+v.^2).^0.5.*(5)); %向上延拓5米
%频率域处理
Sp=Sgex(1:nyn2,:).*f;
Sp=[Sp;conj(Sp((floor((nyn-1)/2)+1):-1:2,1)) conj(Sp((floor((nyn-1)/2)+1):-1:2,end:-1:2))];
%反变换
gg=ifft2(Sp);
gg=real(gg(eyn+1:end-eyn,exn+1:end-exn));
[x,y]=meshgrid(0:5:290,0:5:350);
figure();
contourf(x,y,g+bk);
colorbar;
axis equal;
xlabel('East-west direction');
ylabel('North-south direction');
figure();
% subplot(1,2,1);
contourf(x,y,gg);
colorbar;
axis equal;
xlabel('East-west direction');
ylabel('North-south direction');
% f1=exp(2*pi*(u.^2+v.^2).^0.5.*(-50)); %向上延拓5米
% %频率域处理
% Sp1=Sgex(1:nyn2,:).*f1;
% Sp1=[Sp1;conj(Sp1((floor((nyn-1)/2)+1):-1:2,1)) conj(Sp1((floor((nyn-1)/2)+1):-1:2,end:-1:2))];
%
% %反变换
% gg1=ifft2(Sp1);
% gg1=real(gg1(eyn+1:end-eyn,exn+1:end-exn));
% subplot(1,2,2);
% contourf(x,y,gg1);
% colorbar;
% axis equal;
% xlabel('East-west direction');
% ylabel('North-south direction');
评论2