clear
clc
% ==========================================================
% Gprmax函数用于读取GprMax2D与GprMax3D软件仿真
% 探地雷达模型生成的二进制波形数据
% ==========================================================
% ==========================================================
% ******************有目标的数据****************
% ==========================================================
% % GprMax仿真的目标和墙体的场景图形%%%%%
tic;
filegeo = 'orin.geo';
[meshdata,header,media] = gprmax2g(filegeo);
figure(1);
[MM,NN]= size(meshdata);
%mesh: 存储模型数据，MxN的数组，其中M为Y轴方向的Yee单元数目,N为X轴方向的Yee单元数目；
%M=header.nx=7/0.01=700，header.nx:模型在X轴的偏移次数；N=header.ny=5/0.01=500; header.ny: 模型在 Y 轴的偏移次数；
imagesc((1:NN)*header.dx,(1:MM)*header.dy,meshdata);%header.dx: 模型在 X 轴每次偏移大小(单位:m);
%header.dy: 模型在 Y 轴每次偏移大小(单位:m)
colormap('jet')
axis('xy')
set(gca,'FontSize',18,'FontWeight','Bold')
PT = xlabel('方位向(m)'); set(PT,'FontSize',18,'FontWeight','Bold');
PT = ylabel('距离向(m)'); set(PT,'FontSize',18,'FontWeight','Bold');
PT = title('目标模型'); set(PT,'FontSize',12);
% % GprMax仿真输出目标和墙体的总回波信号%%%%
%Header：该变量包括以下成员：Header.title: 模型的名称;迭代次数Header.iterations=ceil(Header.removed/Header.dt)朝向正无穷取整 ；
%Header.removed=40e-9；(等于#tx:命令的第五个参数，时间窗)
%Header.dx: 在 X 轴每次偏移大小；Header.dy: 在 Y 轴每次偏移大小；
%Header.dt=1/(c*sqrt(1/Header.dx^2+1/Header.dy^2))；最大允许时间步长
%Header.NSteps: 仿真次数；等于*.in 文件中#analysis:命令的第一个参数，天线个数；
fileout = 'orin.out';
[Header,Fields]=gprmax(fileout);
N = 1:Header.NSteps; %仿真次数=天线个数=移动次数，Header.TxStepX=0.01/Header.dx
Position = Header.dx * Header.tx+(N-1)*(Header.dx*Header.TxStepX); %天线每次所在位置%Header.tx=0.85/Header.dx=85；
%Header.ty=2.5/Header.dy=250；等差数列的通项公式
Data_with_target(:,:) = Fields.ez(:,1,:); %转换数组格式
t = Fields.t ;%Fields.t: 每个波形的时间轴。数组大小 Header.iterations*1，迭代次数*1
dt = Header.dt;%最大允许时间步长
% ==========================================================
% ***********消除背景的影响,得到目标的回波信号**********
% ==========================================================
Data_with_target(1:300,:) = 0;%268
Data_target = Data_with_target';
%==========================================================
% 墙体参数 厚度d 相对介电常数 e 以及 电磁波传播速度 c
%==========================================================
d_wall_Low = 0.6;
d_wall_Upper = 0.8;
d = d_wall_Upper - d_wall_Low;
er_concret = 4.5;
speed = 3e8;
%
% %----收发共置天线------------
Number_trace = Header.NSteps;
recx = Header.dx * Header.tx+(N-1)*(Header.dx*Header.TxStepX); %发射天线x位置Header.tx=0.85/Header.dx=85,Header.TxStepX=0.0/Header.dx
recz = Header.dy * Header.ty+(N-1)*(Header.dy*Header.TxStepY);%发射天线y位置Header.ty=2.5/Header.dy=250,Header.TxStepY=0.1/Header.dy
dt1 = Header.dt;%Header.dt=1/(c*sqrt(1/Header.dx^2+1/Header.dy^2))
%***************引入人为噪声****************************
%*******************************************************
Rx_size_array = Number_trace;%= Header.NSteps为20个
[zz,Number_point_signal]=size(Data_target);%注意，这里矩阵转置过了前者是天线个数20，后者是频点数1272
Rec_wave = zeros(Rx_size_array,Number_point_signal);%20*1272
xxxx = 0;
%将Data_target中所有元素平方后加起来，除以频点的个数
for ii = 1 : Rx_size_array %1：20 20个天线
xxxx = xxxx + sum(Data_target(ii,:) .^ 2) / Number_point_signal; %Data_target：20*1272
end
E_Tw = xxxx / Rx_size_array; % 平均信号能量=将Data_target中所有元素平方后加起来，除以频点的个数，再除以天线数（除以二者实际上就是矩阵元素总数）
%这句话的意思就是将Data_target中所有元素平方后加起来除以Data_target中元素总数
SNR = 15; %SNR setup
EbNo = 10.^( SNR ./10); %31.6228
N0 = E_Tw ./ EbNo; % 计算噪声功率
sigma = sqrt(N0./2 ); %计算标准差
for ii = 1 : Rx_size_array %1:24
%产生一个随机分布的指定均值和方差的矩阵：将randn产生的结果乘以标准差，然后加上期望均值即可。
mm = sigma .* randn(1,Number_point_signal); %randn是产生标准正态分布的伪随机数，这句话是产生均值为0，方差为N0./2 的正态分布的伪随机数
Rec_wave(ii,:) = Data_target(ii,:) + mm; %带噪信号 产生的噪声：均值为0，方差为N0./2
% Rec_wave(ii,:) = Data_target(ii,:); %无噪声信号
end
% =================================================================================
% 产生gausian(高斯)脉冲信号
% =================================================================================
% GprMax 的高斯脉冲信号表达式 I = exp(-b(t-a)^2),b = 2 * pi * pi * f^2 , a = 1/f; 1GHZ 幅度为1V(参考user guide)
f = 1e9;
b = 2 * pi^2 * f^2 ;
a = 1 / f;
t = 0 :dt1 : 30e-9;%要注意每次如果重新设置了，这个时间窗一定要改
pulse = exp(-b*(t-a).^2);
% pulse0 = exp(-b*(t-a).^2);
% bipulse = diff(pulse0);
% p1=size(pulse0,1);
% pulse1=zeros(1,p1);
% pulse=[bipulse 0]*10;
%===================================================
%--------成像的参数----------------------
%===================================================
grid_resolution = 0.05; %将目标成像空间离散化
Image_Region_X_Start = 0; Image_Region_X_end = 4;%在X方向上的距离范围
Image_Region_Y_Start = 0.5; Image_Region_Y_end = 2.5;%在Y方向上的距离范围
Image_Region_X = [Image_Region_X_Start : grid_resolution : Image_Region_X_end];
Image_Region_Y = [Image_Region_Y_Start : grid_resolution : Image_Region_Y_end];
Dim_Image_X = length(Image_Region_X);%41
Dim_Image_Y = length(Image_Region_Y);%41
Dim_Image = Dim_Image_X * Dim_Image_Y ;%整个场景划分的点数
% % % 定义降维后的观测数据的性质 Q1=200，Q2=20，,M=1272，N=20 Nx=41 Nz=41
Number_measure =200; %控制采样点数 所谓的Q1
% 产生随机孔径 %选择一部分天线用来成像，其目的还是用来降低使用数据的量，减少后端数据量的压力
rand1 = randperm(Rx_size_array);%产生20维的随机排列的行向量，元素是1:20,原始的天线数20
Number_Aperture=21;%随机选择
index = rand1(1 : Number_Aperture);%
% Number_Aperture=3;%有目的性选
% index =(10:1:12);
Random_apeture = sort(index);%默认把向量元素从小到大排序
Pusi = zeros(Number_Aperture * Number_point_signal,Dim_Image); %dictionary matrix
Phi = zeros(Number_Aperture * Number_measure,Number_Aperture* Number_point_signal);
t1 = zeros(Number_Aperture * Number_measure,1);%Q1Q2*1=3000*1
distance=1;
kk=0;
%-----------开始构造大矩阵---------------
tic;
for mm = 1 : Number_Aperture %1---15
ii = Random_apeture(mm);%确定随机孔径点位置 1 2 3 4 6 8 11 12 13 16 18 19 21 22 23
annte_loca = recx(:,mm); % recx是一个行向量，1.5:0.1:3.4，这一步是在后面正对天线墙体被选择，侧视的墙体被屏蔽
for pos_y = 1: Dim_Image_Y %1:41
Image_Y = Image_Region_Y(pos_y);
for pos_x = 1: Dim_Image_X %1：41
Image_X = Image_Region_X(pos_x);
% 原始字典 0
% theta = atan(abs(recx(ii)-Image_X)/abs(recz(ii)-Image_Y));
% Image_tao = 2 * ((sqrt((recx(ii) - Image_X).^2 + (recz(ii) -Image_Y).^2)) / speed+(d*sqrt(er_concret-sin(theta).^2)-d*cos(theta))/speed);%聚焦时延
% %% ***********************************************************************************************************************************************
% 原始字典
% Image_tao = 2 * (sqrt((recx(ii) - Image_X).^2 + (recz(ii) - Image_Y).^2)) / speed;
% %*********************************************************************************************************************************************
%%%% 原始字典+
Image_tao = 2 * (sqrt((recx(ii) - Image_X).^2 + (recz(ii) - Image_Y).^2)) / speed+2*(d*sqrt(er_concret)-d)/speed;%聚焦时延
% %*****************************************************
