Matlab雷达回波数据模拟.doc

clear, hold offformat compactJ = sqrt(-1);close all% Get root file name for saving resultsfile=input(Enter root file name for data and listing files: ,s);% form radar chirp pulseT = 10e-6; % pulse length, secondsW = 10e6; % chirp bandwidth, Hzfs = 12e6; % chirp sampling rate, Hz; oversample by a littlefprintf(nPulse length = %g microsecondsn,T/1e-6)fprintf(Chirp bandwidth = %g Mhzn,W/1e6)fprintf(Sampling rate = %g Msamples/secn,fs/1e6)s = git_chirp(T,W,fs/W); % 120-by-1 arrayplot(1e6/fs)*(0:length(s)-1),real(s) imag(s)title(Real and Imaginary Parts of Chirp Pulse)xlabel(time (usec)ylabel(amplitude)gridNp = 20; % 20 pulsesjkl = 0:(Np-1); % pulse index array, 慢时间采样的序列，注意第一个PRI标记为0是为了慢时间起始时刻从零开始PRF = 10.0e3; % PRF in HzPRI = (1/PRF); % PRI in secT_0 = PRI*jkl; % relative start times of pulses, in secg = ones(1,Np); % gains of pulsesT_out = 12 40*1e-6; % start and end times of range window in sec， 这个就是接收窗的时间宽度TrecT_ref = 0; % system reference time in usec， T_ref = 0指T_0=0时，r_at_T_0 = ri ；当T_0 = 0时，r_at_T_0 = ri - vi*T_0(j)fc = 10e9; % RF frequency in Hz; 10 GHz is X-bandfprintf(nWe are simulating %g pulses at an RF of %g GHz,Np,fc/1e9)fprintf(nand a PRF of %g kHz, giving a PRI of %g usec.,PRF/1e3,PRI/1e-6)fprintf(nThe range window limits are %g to %g usec.n, . T_out(1)/1e-6,T_out(2)/1e-6)% Compute unambiguous Doppler interval in m/sec% Compute unambiguous range interval in metersvua = 3e8*PRF/(2*fc); %第一盲速rmin = 3e8*T_out(1)/2;rmax = 3e8*T_out(2)/2;rua = 3e8/2/PRF;fprintf(nThe unambiguous velocity interval is %g m/s.,vua)fprintf(nThe range window starts at %g km.,rmin/1e3)fprintf(nThe range window ends at %g km.,rmax/1e3)fprintf(nThe unambiguous range interval is %g km.nn,rua/1e3)% Define number of targets, then range, SNR, and% radial velocity of each. The SNR will be the actual SNR of the target in% the final data; it will not be altered by relative range.Ntargets = 4;del_R = (3e8/2)*( 1/fs )/1e3; % in kmranges = 2 3.8 4.4 4.4*1e3; % in kmSNR = -3 5 10 7; % dBvels = -0.4 -0.2 0.2 0.4*vua; % in m/sec% From SNR, we compute relative RCS using the idea that SNR is proportional% to RCS/R4. Students will be asked to deduce relative RCS.rel_RCS = (10.(SNR/10).*(ranges.4);rel_RCS = db(rel_RCS/max(rel_RCS),power)fprintf(nThere are %g targets with the following parameters:,Ntargets)for i = 1:Ntargets fprintf(n range=%5.2g km, SNR=%7.3g dB, rel_RCS=%7.3g dB, vel=%9.4g m/s, . ranges(i)/1e3,SNR(i),rel_RCS(i),vels(i) )end% Now form the range bin - pulse number data mapdisp( )disp( )disp(. forming signal component)y = radar(s,fs,T_0,g,T_out,T_ref,fc,ranges,SNR,vels); % y是337-by-20的矩阵% add thermal noise with unit powerdisp(. adding noise)%randn(seed,77348911);My,Ny = size(y);nzz = (1/sqrt(2)*(randn(My,Ny) + J*randn(My,Ny); % 产生方差为1的复高斯白噪声y = y + nzz;%create log-normal (ground) clutter with specified C/Nand 具体原理不清楚，需要时套用此格式即可！% log-normal standard deviation for amplitude, uniform phase% Clutter is uncorrelated in range, fully correlated in pulse #disp(. creating clutter)CN = 20; % clutter-to-noise ratio in first bin (dB)SDxdB = 3; % in dB (this is NOT the sigma of the complete clutter)ncc=10 .(SDxdB*randn(My,Ny)/10);ncc = ncc.*exp( J*2*pi*rand(My,Ny) );% Force the power spectrum shape to be Gaussiandisp(. correlating and adding clutter)G = exp(-(0:4).2/1.0);G = G;zeros(Ny-2*length(G)+1,1);G(length(G):-1:2);for i=1:My ncc(i,:)=ifft(G.*fft(ncc(i,:);end% rescale clutter to have desired C/N ratiopcc = var(ncc(:);ncc = sqrt(10(CN/10)/pcc)*ncc;% 10*log10(var(ncc(:)/var(nzz(:) % check actual C/N% Now weight the clutter power in range for assume R2 (beam-limited) losscweight = T_out(1)*(T_out(1) + (0:My-1)*(1/fs).(-1);cweight = cweight*ones(1,Np);ncc = ncc.*cweight; % var(ncc)可以看出20列clutter的方差均在30左右y = y + ncc;My,Ny=size(y);d=(3e8/2)*(0:My-1)*(1/fs) + T_out(1)/1e3; % T_out(1)/1e3是接收窗的起始时刻plot(d,db(y,voltage)xlabel(distance (km)ylabel(amplitude (dB)grid% Save the data matrix in specified file.% Save the student version in the mystery file.% Also save all parameter value displays in corresponding filedata_file=file,.mat;mystery_file=file,_mys.mat;listing_file=file,.lis;eval(save ,data_file, J T W fs s Np PRF PRI T_out fc vua, . rmin rmax rua Ntargets ranges vels SNR rel_RCS y);eval(save -v6 ,mystery_file, J T W fs s Np PRF T_out fc y);fid=fopen(listing_file,w);fprintf(fid,rDESCRIPTION OF DATA IN FILE ,file,.mat AND ,file,_mys.matrr);fprintf(fid,rPulse length = %g microsecondsr,T/1e-6);fprintf(fid,Chirp bandwidth = %g Mhzr,W/1e6);fprintf(fid,Sampling rate = %g Msamples/secr,fs/1e6);fprintf(fid,rWe are simulating %g pulses at an RF of %g GHz,Np,fc/1e9);fprintf(fid,rand a PRF of %g kHz, giving a PRI of %g usec.,PRF/1e3,PRI/1e-6);fprintf(fid,rThe range window limits are %g to %g usec.r, . T_out(1)/1e-6,T_out(2)/1e-6);fprintf(fid,rThe unambiguous velocity interval is %g m/s.,vua);fprintf(fid,rThe range window starts at %g km.,rmin/1e3);fprintf(fid,rThe range window ends at %g km.,rmax/1e3);fprintf(fid,rThe unambiguous range interval is %g km.rr,rua/1e3);fprintf(fid,rThere are %g targets with the following parameters:, . Ntargets);for i = 1:Ntargets fprintf(fid,r range=%5.2g km, SNR=%7.3g dB, rel_RCS=%7.3g dB, vel=%9.4g m/s, . ranges(i)/1e3,SNR(i),rel_RCS(i),vels(i) );endfclose(fid);fprintf(nnData is in file ,data_file)fprintf(nStudent data is in file ,mystery_file)fprintf(nListing is in file ,listing_file,nn)_用到的函数function y = radar( x, fs, T_0, g, T_out, T_ref, fc, r, snr, v )% RADAR simulate radar returns from a single pulse or burst% of identical pulses% usage:% R = radar( X, Fs, T_0, G, T_out, T_ref, Fc, R, SNR, V )% X: baseband single pulse waveform (complex vector)% Fs: sampling frequency of input pulse in Hz% T_0: start time(s) of input pulse(s) sec% (number of pulses in burst assumed = length(g) )% G: complex gain(s) of pulse(s)， 即慢时间，各个PRI对应的脉冲的前的加权 20-by-1% T_out: 2-vector T_min,T_max defines output% window delay times w.r.t. start of pulse% T_ref: system reference time, needed to simulate% burst returns. THIS IS THE t=0 TIME !% Fc: center freq. of the radar. in Hz% R: vector of ranges to target(s) meters% (number of targets assumed = length(r) )% SNR: vector of target SNRs (unit noise power assumed)% This will be SNR *after* allowing for R4% V: vector of target velocities (optional) in m/sec% (positive velocities are towards the radar)% note(1): VELOCITY in meters/sec !% distances in m, times in sec, BW in Hz.% note(2): assumes each pulse is constant (complex) amplitude% note(3): will accomodate up to quadratic phase pulses% note(4): vector of ranges, R, allows DISTRIBUTED targets% (c) jMcClellan 7/28/90% Modified by M. A. Richards, August 1991J = sqrt(-1);c = 3e8; % velocity of light in m/secMx = length(x);delta_t = 1/fs; % sampling interval (sec)t_y = T_out(1):delta_t:T_out(2) ; % output sampling times (sec)，接收窗的宽度内的等间隔采样 337-by-1T_p = Mx*delta_t; % length of input pulse (sec)，基带信号chirp的脉冲持续时间，即Te% Assume zero velocities (stationary targets) if no velocity% vector providedif nargin 7 v = zeros(r);end% ensure that all vectors are column vectorsx=x(:); g=g(:); T_0=T_0(:); r=r(:); snr=snr(:); v=v(:);% determine the quadratic phase modulation parameters for% later interpolation of pulse samplest_x = delta_t*0:(Mx-1);x_ph = unwrap(angle(x); %基带chirp信号的相位，可以看出x_ph是个抛物线q = polyfit(t_x,x_ph,2); %目的是用 q = ax2 + bx + c 逼近x_ph% check result using correlation coefficientxfit = polyval(q,t_x); % 看看用 q = ax2 + bx + c 拟合的xfit与x_ph的一致程度if (x_ph*xfit)/norm(x_ph)/norm(xfit) = T_out(2) | tmax =0 & t_vals T_p); %T_p是chirp基带信号的长度，一个chirp脉冲携带有效测量数据，即Te上的采样点 if tau T_out(2) fprintf(nEcho from target #%g at range %g km,i,ri) fprintf(nFINISHES AFTER the range window) fprintf(non pulse#%g.n,j) end % Place scaled, range-delayed, Doppler shifted pulse into output matrix % Unit noise power and unit nominal pulse amplitude assumed to % get amplitude from SNR. amp = 10(snr(i)/20);% n_out 是对应chirp脉冲宽度Te的120-by-1向量，原来接收的chirp信号未经过脉冲压缩在距离上占据c*Te/2米的长度，和对应长度的rect信号在距离上占据的长度是一样的! 经过脉冲压缩后chirp信号在距离上才占据c/(2B)米的长度。 y(n_out,