dopplercanshu.rar_Radon参数估计_radar_radon 估计_radon 多普勒_雷达Radon变换

clear all; clc; %信号载波 c=3.e8; %光速 fc=1.e10; %信号载波频率 lamda=c/fc; %波长 %发射信号参数 Tp=10.e-6; %信号脉宽 Br=100.e6; %信号带宽 kr=Br/Tp; %线性调频斜率 fs=200.e6; %采样频率 pr=c/(2*fs); %距离向像素间隔 Nr=ceil(fs*Tp); %距离向采样点数 Tr=Nr/fs; %雷达运动几何参数 H=6.e3; %雷达高度 V=100.; %雷达飞行速度 Rg=30.e3; %雷达正侧视距离 bias=2/180*pi; %接收斜视角 R0=sqrt(Rg^2+H^2)/cos(bias); %斜距 Ls=320.; %合成孔径长度 Ts=Ls/V; %合成孔径时间 PRF=200.; %脉冲重复频率 Na=ceil(PRF*Ts); %方位向采样点数 Ta=Na/PRF; tr=2*R0/c+linspace(-Tr/2,Tr/2,Nr); %快时间 ta=linspace(-Ta/2,Ta/2,Na); %慢时间 fdc0=2.*V.*sin(bias)/lamda; %多普勒中心频率理论值 fdr0=-2.*V^2/lamda/R0; %多普勒调频率理论值 %雷达位置 xr=R0*sin(bias)-V*ta; yr=-Rg; %仿真点目标位置 xm=[0.]; ym=[0.]; % xm=[7.,7., 7., 7., 7.,5.25,3.5,1.75, 0.,-1.75,-3.5,-5.25,-7.,-7.,-7.,-7.,-7.,-7., -7., -7., -7.,-5.25,-3.5,-1.75, 0.,1.75, 3.5,5.25, 7., 7., 7., 7.,0.]; % ym=[0.,5.,10.,15.,20., 20.,20., 20.,20., 20., 20., 20.,20.,15.,10., 5., 0.,-5.,-10.,-15.,-20., -20.,-20., -20.,-20.,-20.,-20.,-20.,-20.,-15.,-10.,-5.,0.]; figure(1) scatter(ym,xm,'.'); %绘出仿真目标场景图 % axis([-25 25 -10 10]); xlabel('距离向'); ylabel('方位向'); title('仿真点目标位置关系'); % ==========================二维回波信号=============================== n=length(xm); %目标数 data=zeros(Nr,Na); for k=1:Na m=1; record=zeros(1,Nr); while m<=n Rt=sqrt((xm(m)-xr(k))^2+(ym(m)-yr)^2+H^2); record=record+exp(-j*4*pi*Rt/lamda).*exp(j*pi*kr*(tr-2*Rt/c).^2).*(abs(tr-2*Rt/c)<Tp/2); m=m+1; end data(:,k)=record; end figure(2) imagesc(abs(data)); xlabel('方位向采样点数'); ylabel('距离向采样点数'); % figure; % datafft=fft(sum(data),1024); % x=linspace(-PRF/2,PRF/2,1024); % plot(x,fftshift(abs(datafft))); % ============================距离压缩============================== fr=linspace(-fs/2,fs/2,Nr); %距离向频率 fr=fftshift(fr); %fr左右两半对换 ref=exp(j*pi*fr.^2/kr); %距离向滤波器 for k=1:Na data(:,k)=ifft(fft(data(:,k)).*ref.'); end figure(3) imagesc(1:Na,1:Nr,abs(data)); title('距离压缩'); xlabel('方位向'); ylabel('距离向'); %%%%%%%%%%%%%%%%%%%%%%%%%多普勒中心频率估计%%%%%%%%%%%%%%%%%%%%%%%%%%% % ===========================Radon变换============================== a=6; %矩形窗口大小 datax=data(floor(Nr/2)-a:floor(Nr/2)+a,:); %截取待估计数据平面 theta=0:0.19:179; %theta表示Radon平面横轴变量 [amp,rho]=radon(datax,theta); %调用radon函数进行Radon变换 %amp表示Radon变换结果的幅值矩阵 %rho表示Radon平面的纵轴变量 max_a=max(max(abs(amp))); figure(4) mesh(theta,rho,abs(amp)/max_a); xlabel('\theta'); ylabel('\rho'); zlabel('归一化幅度'); grid on; % =========================多普勒中心频率值======================= nn=6; [row,col]=size(amp); %row为amp矩阵的行数,col为列数 ampn=amp.^nn; %待积分数组矩阵 for l=1:col radon1_int(l)=sum(ampn(:,l)); end figure(5) plot(theta,abs(radon1_int/max(radon1_int))); xlabel('\theta'); ylabel('归一化幅度'); [radon1_int_max,index1]=max(radon1_int); theta_est1=theta(index1); %amp_int最大幅值对应的角度 de_r=c/fs; K=PRF*de_r; %变换比例系数 fdc=cot(theta_est1*pi/180)*K/lamda %多普勒中心频率估计值 %%%%%%%%%%%%%%%%%%%%%%%%%%多普勒调频率估计%%%%%%%%%%%%%%%%%%%%%%%%%%% % ==========================距离走动校正============================= deta_r=lamda*fdc*ta; %距离走动 datam=zeros(Nr,Na); for k=1:Na datam(:,k)=ifft(fft(data(:,k)).*exp(j*2*pi*fr*deta_r(k)/c)'); end figure(6); imagesc(abs(datam)); %距离走动校正后图像 title('距离走动校正'); % ===========================Wigner变换============================== data1=datam(500,:); %选取一个距离单元的数据 t0=linspace(1,Na,Na); figure(7) plot(t0,real(data1)); %该距离单元信号波形 [result,t1,f1]=tfrspwv(data1'); %调用tfrspwv工具箱函数 result=real(result); figure(8) imagesc(f1,t1,result); %时频分析图 % ===========================radon变换=============================== theta=0:1:179; [result1,xp]=radon(result,theta); result1_max=max(max(result1)); %最大值 result1=result1/result1_max; %幅值归一化 figure(9) mesh(theta,xp,abs(result1)); % ===========================多普勒调频率值========================= [m,n]=size(result1); resultn=result1.^nn; for ii=1:n radon2_int(ii)=sum(resultn(:,ii)); end figure(10) plot(theta,abs(radon2_int/max(radon2_int))); xlabel('\theta'); ylabel('归一化幅度'); [radon2_int_max,index2]=max(radon2_int); theta_est2=theta(index2); %积分曲线最大幅值对应的角度 fdr=-cot(theta_est2*pi/180)/lamda %调频率估计值