mytp.m

这是对波场做的一个TP变换

load('z.dat');
[M,N]=size(z);
for i=1:N
 uxt(1:M,i)=z(1:M,i);
end
poinum=M;
tranum=N;
poinum2=M;
tranum2=N;
deltx=10.0;
dx=1.0;
dt=0.0001;
pmax=0.01;
pmin=0.0;
dp=0.00005;
pnum=round((pmax-pmin)/dp)+1;
utp(1:poinum,1:pnum)=0;
for luoi=1:poinum 
    for luoj=1:pnum   
        for luok=1:tranum 
            offset=deltx+(luok-1)*dx;  
            tt=((luoj-1)*dp+pmin)*offset/dt;
            tt=round(luoi+tt);
            if tt>0 && tt<=poinum
                utp(luoi,luoj)=utp(luoi,luoj)+uxt(tt,luok);
            end
        end
    end
end

uxt2(1:poinum2,1:tranum2)=0;
for luoi=1:poinum2
    for luoj=1:tranum2
        offset=deltx+(luoj-1)*dx;
        for luok=1:pnum
            tt=((luok-1)*dp+pmin)*offset/dt;
            tt=round(luoi-tt);
            if tt>0 && tt<=poinum2
                uxt2(luoi,luoj)=uxt2(luoi,luoj)+utp(tt,luok);
            end
        end
    end
end
for luoi=1:tranum2
    maxvalue=max(uxt2(:,luoi));
    for luoj=1:poinum2
        uxt2(luoj,luoi)=uxt2(luoj,luoi)/maxvalue;  %uxt2最后输出的很小，和输入的z相差很大，相对值差不多
    end
end
save ztp.dat utp -ascii;   %保存TP图
save z2.dat uxt2 -ascii;   %保存反变换回去的图
% for i=1:N
%  xt11(1:M,i)=z(1:M,i)-uxt2(1:M,i);
% end
% save z11.dat xt11 -ascii;  %保存输入数据和反变换回去的数据差值图