hcgpr.m

雷达正演程序

function hcgpr
clear;
prompt={'电性分界面的个数'}; title='合成雷达记录的计算';
defaults={'3 '};
cells=inputdlg(prompt,title,1,defaults);
if isempty(cells)
  warndlg('您选择了对话框中的“Cancel" 键，退出程序运行',title);
  return;
end
n=str2num(cells{1});
prompt={'第一层的厚度','第二层的厚度','第三层的厚度','第四层的厚度'};
defaults={'50 ','30','50','30'};
cells=inputdlg(prompt,title,1,defaults);
if isempty(cells)
  warndlg ('您选择T对话框中的“Cancel"键，退出程序运行',title);
  return;
end
for i=1:n
h(i)=str2num(cells{i});
end
prompt={'第一层的相对介电常数','第二层的相对介电常数','第三层的相对介电常数','第四层的相对介电常数','第五层的相对介电常数'};
defaults={'1' ,'6','9','12','20'};
cells=inputdlg(prompt,title,1, defaults);
if isempty(cells)
warndlg ('您选择了对话框中的“Cancel”键，退出程序运行',title);
return;
end
for i=1:(n+1)
E(i)=str2num(cells{i});
end
prompt={'第一层的电阻率','第二层的电阻率','第三层的电阻率','第四层的电阻率'};
defaults={'100000000','5000','2500','2000'};
cells=inputdlg(prompt,title,1, defaults);
if isempty(cells)
warndlg ('您选择了对话框中的“Cancel”键，退出程序运行',title);
return;
end
for i=1:n
p(i)=str2num(cells{i});
end
c=30;
for i=1:n
v(i)=c./sqrt(E(i));
end
B(1)=0;
for i=1:n
   if p(i)==100000000  B( i)=0;
   else B(i)=188*(1/p(i))/sqrt(E(i));
   end
end
for j=1:n
A(j)=exp(-2*B(j).*h(j)./100);
end
for i=1:n
R(i)=[sqrt(E(i))-sqrt(E(i+1))]/[sqrt(E(i))+sqrt(E(i+1))];
end
r(1)=R(1); a=1; b=1; c=0;
for i=2:n
a=a*A(i);b=b*[1-R(i-1)^2];c=c+h(i);r(i)=a*R(i)*b/[2*c/100];
end
for i=1:n
t(i)=0;
  for j=1:i
      x(j)=0; x(j)=2*h(j)./v(j);t(i)=t(i)+x(j);
  end
end
prompt={'雷达天线的中心频率(MHz)','所计算的时间坐标长度(ns)','理论计算时的时间点距(ns)'}
defaults={'1000','30','0.01'};
cells=inputdlg(prompt,title,1,defaults);
if isempty(cells)
warndlg('您选择7对话框中的“Cancel" 键，程序退出运行',title)
return;
end
fl=str2num(cells{1}); f=fl./1000; tt=str2num(cells{2});
dt=str2num(cells{3});
ii=1;
for T=0:dt:tt
z1=sin(2*pi*T*f);z2=exp(2*pi*T*f/sqrt(3)); z3=T.*T;
y(ii)=z3.*z1./z2; ii =ii+1;
end
for i=1:n
t(i)=round(t(i).*(1/dt))*dt;
end
for k=1:(tt/dt+1)
RR(k)=0;
end
for i=1:n
    j(i)=round(t(i).*(1/dt)+1);RR(j(i))=r(i);
end
k=1:(tt/dt+1);RR(k ); c=conv(RR(k),y(k ));
data_store=questdlg('您想存储合成雷达数据结果吗?',title);
if strcmp(data_store,'YES')
cells=inputdlg('存储合成雷达结果',title,1,{'gpr.txt'});
if isempty(cells)
warndlg ('您选择了对话框中的“Cancel"键，程序退出运行',title);
return;
end
file= cells{1};fid=fopen(file,'w ' ); k=1:(tt/dt+1);
z=[(k -1)*dt;c(k)];fprintf(fid,'%6.2f% g\t',z);
status= fclose(fid);
else
pause(1);
end
inter=max(c(k)); prompt={'请输入数据道数目'};defaults={'30'};
dells=inputdlg(prompt,title,1, defaults);traces=str2num(dells{1});
stype=menu('选择所要绘制的曲线类型','绘制合成雷达记录Swiggle曲线','绘制合成雷达记录平面等值线');
prompt={'第一个界面反射信号的增益','第二个界面反射信号的增益','第三个界面反射信号的增益','第四个界面反射信号的增益'};
defaults={'1.0','1.0','1.0','1.0'};
fell=inputdlg(prompt,title,1, defaults);
if isempty(fell)
warndlg ('您选择T对话框中的“Cancel”键，程序退出运行',title);
return;
end
for i=1:n
apply(i)=str2num(fell{i});
end
for i=1:n-1
t_inter(i)=round((t(i)+(t(i+1)-t(i))./2)./dt);
end
i=2:n-1;
for k=1:(tt/dt+1)
   if k<=t_inter(1); d(k)=c(k).*apply(1);
   elseif k>t_inter(i); d(k)=c(k).*apply(i);
   end
end
k=1:(tt/dt+1);
for i=1:traces
dd(i,k)=d(k);
end
switch stype
case 1
figure(1); hdzqx(tt,dt,traces,inter,dd);
case 2
figure(1); hdzx(traces,tt,dt,dd);
end
for i=1:traces
   cc(i,k)=c(k);
end
promp={'输入噪声系数'};defaults={'0.25'};
gells=inputdlg(promp,title,1,defaults);zs=str2num(gells{1});
Max=log(max(c(k)));M=fix(Max-1);
sj=randn(traces,tt/dt+1).*(10-M).*zs;
for i=1:traces
  for k=1:tt/dt+1
cc_sj(i,k)=cc(i,k)+sj(i,k);
end
end
stype=menu('选择所要绘制的曲线类型',' 加噪声后的合成雷达记录Swiggle63曲线','加噪声后的合成雷达记录平面等值线');
switch stype
case 1
figure(2);  hdzqx(tt,dt,traces,inter,cc_sj);
case 2
figure(2);  hdzqx(traces,tt,dt,cc_sj);
end
prompp={'傅立叶变换数据量参数'};defAns={'10'};
hells=inputdlg(prompp,title,1,defAns); tn=str2num(hells{1});
filter_Freq(tn,tt,dt,traces,inter,cc_sj );
function hdzqx(tt,dt,traces,inter,aa )
k=1:tt/dt+1; 
hold on;
for i=1:traces
plot(aa(i,k)./inter+i-1,(k-1)*dt);
end
axis([-1/2 traces-1/2 0 tt]); axis('ij');
xlabel('道号');ylabel('双程旅行时 t(ns)');
hold off;
function hdzx(traces,tt,dt,aa)
i=1:traces; k=1:(tt/dt+1);z=aa(i,k);
contour(i-1,(k-1).*dt,z); axis('ij');
xlabel('道号');ylabel('双程旅行时 t(ns)');
function filter_Freq(tn,tt,dt,traces,inter,cc_sj)
k1=tt/dt+1; k2=1:traces;transf=fft2(cc_sj);
filttransf(k2,1:tn)=transf(k2,1:tn);
filttransf(k2,(k1-tn):k1)=transf(k2,(k1-tn):k1);
filtsig=ifft2(filttransf); signal=real(filtsig);
type=menu('选择所要绘制的曲线类型','滤波后的Swiggle曲线','滤波后的平面等值线');
switch type
case 1
figure(3); hdzqx(tt,dt,traces,inter,signal);
case 2
figure(3); hdzx(traces,tt,dt,signal);
end