radon_radon变换目标_Radon_Radon变换_

clc clear all close all %%%%%%%%%%%%%%%%%%设定运动目标运动参数%%%%%%%%%%%%%%%%%%%%%%%%%%%% v=5; %运动目标的速度 t0=30; %最近距离对应时刻 c=1500; %水中声速 r0=5000; %最近距离 a=v/(r0*10); theta0=30; %运动目标初始方位角 fs=10000; t=1:60; fL=500; fH=2000; n=fs; fl=round(fL*n/fs)+1; fh=round(fH*n/fs)+1; Wn=2*[fL fH]/fs; bb=fir1(2^12,Wn,'bandpass'); left=1; %1接收点左侧通过，0接收点右侧通过 if left==1 theta_Z =270; %方位解算初始角度设置 theta_k=-1; %用于生成和解算方位的系数 end if left==0 theta_Z =90; %方位解算初始角度设置 theta_k=1; end X0=r0*cosd(theta0+theta_Z); Y0=r0*sind(theta0+theta_Z); x=X0+(t-t0)*v*60*cosd(theta0); y=Y0+(t-t0)*v*60*sind(theta0); theta_XY=mod(round(90-(atan2(x,y)*180/pi)),360); figure plot(t,theta_XY) ylim([0 360]) xlabel('时间/s') ylabel('方位角/度') grid on Dzhixian1=0; Dr=sqrt(Dzhixian1^2+(r0/1000)^2); TL=60+15*log10(Dr); Daisnr=[3 6 9 12 15 18 21 24 ];%最近距离处带内信噪比 SL=Daisnr+TL+70;%船的声源级 figure plot(Daisnr,SL,'rp'); hold on; plot(Daisnr,SL, 'k--.') ; xlabel('带内信噪比/dB'); ylabel('目标声源级/dB'); grid on tt=(1:fs)/fs; snr=zeros(length(Daisnr),length(t)); for mmmm=1:length(Daisnr) for tttt=1:length(t) Dzhixian1=v*60*(t(tttt)-t0)/1000;%Dzhixian1=v*60*(t(tttt)-t0)/1000 Dr1=sqrt(Dzhixian1^2+(r0/1000)^2);%%%%需要人为调距离 SNR1(mmmm,tttt)=SL(mmmm)-70-60-15*log10(Dr1);%%真实的带内信噪比 snr(mmmm,tttt)=10^(SNR1(mmmm,tttt)/20); end end theta=zeros(1,length(t)); %%%%%%%%%%%%%%%%%%%循环%%%%%%%%%%%%%%%%%%%%%%% for mmmm=1:length(Daisnr) mmmm for number=1:100 number theta_scan=0:359; sig=normrnd(0,1,1,length(tt)); sig_narrow=filter(bb,1,sig); s1= sig_narrow/std(sig_narrow);%归一化信号能量 Theta_time=zeros(length(t),length(theta_scan)); %为时间方位矩阵预分配存储空间 for i=1:length(t) t(i)=i; theta(1,i)=mod(-atand(theta_k*(v/r0)*60*(t(i)-t0))+theta_Z+theta0,360); %某时刻下匀速直线运动目标方位 theta(1,i)=mod(-atand(theta_k*(v/r1(mmmm))*60*(t(i)-t0))+theta_Z+theta0,360); end for i=1:length(t) noise_ran1=normrnd(0,1,[1,length(tt)]); %产生3路噪声信号 noise_ran2=normrnd(0,1,[1,length(tt)]); noise_ran3=normrnd(0,1,[1,length(tt)]); noise_narrow1=filter(bb,1,noise_ran1); noise_narrow2=filter(bb,1,noise_ran2); noise_narrow3=filter(bb,1,noise_ran3); noise1= noise_narrow1/std( noise_narrow1); noise2= noise_narrow2/std( noise_narrow2); noise3= noise_narrow3/std( noise_narrow3); Vx(mmmm,:)=(s1*snr(mmmm,i))*cosd(theta(1,i))+1*noise1; %某时刻下运动目标振速 Vy(mmmm,:)=(s1*snr(mmmm,i))*sind(theta(1,i))+1*noise2; %某时刻下运动目标振速 P(mmmm,:)=s1*snr(mmmm,i)+1*noise3; %某时刻下运动目标声压 Vx_fft(mmmm,:)=fft(Vx(mmmm,:),n); Vy_fft(mmmm,:)=fft(Vy(mmmm,:),n); P_fft(mmmm,:)=fft(P(mmmm,:),n); Ix(mmmm,:)=real(conj(Vx_fft(mmmm,:)).*P_fft(mmmm,:)); Iy(mmmm,:)=real(conj(Vy_fft(mmmm,:)).*P_fft(mmmm,:)); for j=fl:fh %%1:length(Vx_fft) theta_est(mmmm,j)=mod(round(atan2(Iy(mmmm,j),Ix(mmmm,j))*180/pi),360); %计算某时刻下每个频点对应的角度 Theta_time(i,theta_est(mmmm,j)+1)=Theta_time(i,theta_est(mmmm,j)+1)+1; end end % figure % pcolor(theta_scan,t,Theta_time) % xlabel('方位角/度') % ylabel('时间/s') % % title('t-theta扫描图'); % colorbar % colormap('jet') % shading interp k_scan=0:a:20*a;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%人为调 theta0_scan=0:1:359; t0_scan=1:length(t); Theta0_k1=zeros(length(t0_scan),length(k_scan),length(theta0_scan)); for k=1:length(t) t00(k)=k; for j=1:length(k_scan) for i=1:length(theta0_scan) for kk=1:length(t) thetaKK=theta_Z-atand(theta_k*60*(t(kk)-t00(k))*k_scan(j))+theta0_scan(i);%% thetaKK=theta_Z-atand(theta_k*60*(t(kk)-t00(k))*k_scan(j))+theta0_scan(i) temp_m11=mod(round(round(thetaKK*10)/10),360)+1; vector11(kk)=Theta_time(kk,temp_m11); end Theta0_k1(k,j,i)=sum(vector11); end end end Theta0_k12=zeros(length(t0_scan),length(k_scan)); for k=1:length(t) for j=1:length(k_scan) Q1=max(Theta0_k1(k,j,:)); Theta0_k12(k,j)=Q1; end end [A,T]=max(max(Theta0_k12')); t0_est=T; Each_t0_est(mmmm,number)=t0_est; % figure % pcolor(k_scan,t0_scan,Theta0_k12) % xlabel('k') % ylabel('t/s') % % title('t0-K扫描图'); % colorbar % colormap('jet') % shading interp Theta0_k13=zeros(length(t0_scan),length(theta0_scan)); for k=1:length(t) for i=1:length(theta0_scan) Q2=max(Theta0_k1(k,:,i)); Theta0_k13(k,i)=Q2; end end [B,E]=max(max(Theta0_k13)); theta0_est=(E-1); Each_theta0_est(mmmm,number)=theta0_est; % figure % pcolor(theta0_scan,t0_scan,Theta0_k13) % xlabel('theta0/度') % ylabel('t/s') % % title('t0-theta0扫描图'); % colorbar % colormap('jet') % shading interp Theta0_k14=zeros(length(k_scan),length(theta0_scan)); for j=1:length(k_scan) for i=1:length(theta0_scan) Q3=max(Theta0_k1(:,j,i)); Theta0_k14(j,i)=Q3; end end [C,K]=max(max(Theta0_k14')); k_est=(K-1)*a; Each_k_est(mmmm,number)=k_est; % figure % pcolor(theta0_scan,k_scan,Theta0_k14) % xlabel('theta0/度') % ylabel('k') % % title('theta0-K扫描图'); % colorbar % colormap('jet') % shading interp %%%%%%%%%%%%%%%%蒙特卡洛循环%%%%%%%%%%%%%%%%%%%%%%% Tt(1,number)=(t0_est-t0)^2; Theta000(1,number)=(theta0_est-theta0)^2; KKK(1,number)=(k_est-v/r0)^2; end T0_est=sqrt(sum(Tt,2)/(number-1)); Theta000est=sqrt(sum(Theta000,2)/(number-1)); KKK0est=sqrt(sum(KKK,2)/(number-1)); FT0_est(1,mmmm)=T0_est; FTheta000est(1,mmmm)=Theta000est; FKKK0est(1,mmmm)=KKK0est; end figure plot(Daisnr,FT0_est,'r*'); hold on; plot(Daisnr,FT0_est, 'k--.','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','LineWidth',0.5) ; xlabel('最近通过距离处信噪比(dB)'); ylabel('t0标准差(秒)'); % set(gca,'xlim',[0 20]); grid on figure plot(Daisnr,FTheta000est,'r*'); hold on; plot(Daisnr,FTheta000est, 'k--.','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','LineWidth',0.5) ; xlabel('最近通过距离处信噪比(dB)'); ylabel('theta0标准差(°)'); % set(gca,'xlim',[0 20]); grid on figure plot(Daisnr,FKKK0est,'r*'); hold on; plot(Daisnr,FKKK0est, 'k--.','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','LineWidth',0.5) ; xlabel('最近通过距离处信噪比(dB)'); ylabel('k标准差(1/s)'); grid on %%%%%%%%%%%%%%%%%%%%%%%测速%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fs=10000; ttT=(1:1*fs*60)/fs; c=1500; t=1:60; t0=30; f0=90; v=5; r0=5000; theta0=30; theta_Z =270; %方位解算初始角度设置 theta_k=-1; %用于生成和解算方位的系数 k=v/r0; fL=40; fH=140; Wn=2*[fL fH]/fs; bb=fir1(2^12,Wn,'bandpass'); for i=1:length(ttT) thetaKKK1(i)=theta_Z-atand(theta_k*(ttT(i)-t0)*60*k)+theta0; end delt_f1=-1*v*60*f0*(cosd(thetaKKK1-theta0))/c; f=(f0+delt_f1);%*30; fl=min(f); fh=max(f); vv=(f-fl)./(fh-fl)*2-1; Sig1=0.5*vco(vv,[fl,fh],fs); s1=Sig1/std(Sig1); Daisnr=[3 6 9 12 15 18 21 24 ];%%最近距离处带内信噪比 TL=60+15*log10(r0/1000); SL=Daisnr+TL+70; for mmmm=1:length(SL) for ttTtt=1:length(ttT) Dzhixian=v*60*(ttT(ttTtt)-t0)/1000; Dr=sqrt(Dzhixian^2+(r0/1000)^2);%%%% SNR1(mmmm,ttTtt)=SL(mmmm)-70-60-15*log10(Dr);%%真实的带内信噪比 snr(mmmm,ttTtt)=10^(SNR1(mmmm,ttTtt)/20); end end for mmmm=1:length(SL) for number=1:100 ttT=1:60; for i