Matlab雷达回波数据模拟

1、clear, hold off format compactJ = sqrt(-1);close all% Get root for saving results(Enter root for data and listing files: ,s);% form radar chirp pulseT = 10e-6;% pulse len gth, sec ondsW = 10e6;% chirp ban dwidth, Hzfs = 12e6;% chirp sampli ng rate, Hz; oversample by a littlefpri ntf(nPulse len gth=%

2、g microsec ondsn ：T/1e-6)fprin tf(Chirp ban dwidth = %g Mhzn,W/1e6)fpri ntf(Sampli ng rate = %g Msamples/sec n,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;% 2

3、0 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 = on es(1,Np);% gains of pulsesT_out = 12 40*1e-6; % start and end times of range window in sec 这个就是接收 窗的时间宽度TrecT_

4、ref = 0;% system referenee 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 inHz; 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)fpri n

5、tf(nThe range wi ndow limits are %g to %g usec.n,.T_out(1)/1e-6,T_out(2)/1e-6)% Compute un ambiguous Doppler in terval in m/sec % Compute un ambiguous range in terval in metersvua = 3e8*PRF/(2*fc); %!盲速 rmin = 3e8*T_out(1)/2;rmax = 3e8*T_out (2)/2;rua = 3e8/2/PRF;fpri ntf(nThe un ambiguous velocity

6、in terval is %g m/s.,vua)fpri ntf(nThe range wi ndow starts at %g km.,rmi n/1e3)fpri ntf(nThe range wi ndow ends at %g km.,rmax/1e3)fprintf(inThe unambiguous range interval is %g km.nn,rua/1e3)% Define number of targets, the n ran ge, SNR, and% radial velocity of each. The SNR will be the actual SNR

7、 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 kmran ges = 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 relativeRCS usi ng the idea that SNR is proportio nal% to R

8、CS/RA4. Stude nts will be asked to deduce relative RCS.rel_RCS = (10A(SNR/10).*(ra nges.A4);rel_RCS = db(rel_RCS/max(rel_RCS),power)fprintf(nThere are %g targets with the following parameters:,Ntargets)for i = 1:Ntargetsfprintf(n range=%5.2g km, SNR=%7.3g dB, rel_RCS=%7.3g dB, vel=%9.4g m/s, .ran ge

9、s(i)/1e3,SNR(i),rel_RCS(i),vels(i)end% Now form the range bin - pulse nu mber data mapdisp()disp()disp(forming sig nal comp onen t)y = radar(s,fs,T_0,g,T_out,T_ref,fc,ranges,SNR,vels); % y是 337-by-20 的矩阵 % add thermal no ise with unit powerdisp(addi ng no ise)%randn( seed,77348911);My,Ny = size(y);n

10、zz = (1/sqrt(2)*(randn(My,Ny) + J*randn(My,Ny); % 产生方差为 1 的复高斯白噪 声y = y + nzz;% create log-normal (ground) clutter with specified C/N and 具体 原理不清楚，需要时套用此格式即可！% log-no rmal sta ndard deviati on for amplitude, uniform phase% Clutter is uncorrelated in range, fully correlated in pulse #disp(creating cl

11、utter)CN = 20;% clutter-to-no ise ratio in first bin (dB)SDxdB = 3;% in dB (this is NOT the sigma of the complete clutter)ncc=10 .A(SDxdB*ra ndn (My,Ny)/10);ncc = n cc.*exp( J*2*pi*ra nd(My,Ny) );% Force the power spectrum shape to be Gaussiandisp(correlating and adding clutter)G = exp(-(0:4).A2/1.0

12、);G = G;zeros(Ny-2*le ngth(G)+1,1);G(le ngth(G):-1:2);for i=1:Myn cc(i,:)=ifft(G.*fft (n cc(i,:);end% rescale clutter to have desired C/N ratiopcc = var(n cc(:);ncc = sqrt(10A(CN/10)/pcc)* ncc;% 10*log10(var(ncc(:)/var(nzz(:) % check actual C/N% Now weight the clutter power in range for assume RA2 (

13、beam-limited) loss cweight = T_out(1)*(T_out + (0:My-1)*(1/fs).A(-1);cweight = cweight* on es(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(dista nee (km)ylabe

14、l(amplitude (dB)grid% Save the data matrix in specified file.% Save the stude nt vers ion in the mystery file.% Also save all parameter value displays in eorresponding 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 fe vua, .rmin rmax ru

15、a Ntargets ranges vels SNR rel_RCS y);eval(save -v6 ,mystery_file, J T W fs s Np PRF T_out fe y);fid=fopen(listing_file,w);fprintf(fid,rDESCRIPTIONOF DATA IN FILE,file,.mat AND ,file,_mys.matr门);fprin tf(fid,rPulse len gth = %g mieroseeo ndsr,T/1e-6);fprin tf(fid,Chirp ban dwidth = %g Mhzr,W/1e6);fp

16、rintf(fid,Sampling rate = %g Msamples/seer,fs/1e6);fprintf(fid,rWe are simulating %g pulses at an RF of %g GHz,Np,fc/1e9);fprin tf(fid,ra nd a PRF of %g kH z, givi ng a PRI of %g usee.,PRF/1e3,PRI/1e-6);fprin tf(fid,rThe range win dow limits are %g to %g usee.r, .T_out(1)/1e-6,T_out(2)/1e-6);fprin t

17、f(fid,rThe un ambiguous veloeity in terval is %g m/s.,vua);fprintf(fid,rThe range window starts at %g km.,rmin/1e3);fprin tf(fid,rThe range window ends at %g km.,rmax/1e3);fprintf(fid,rThe unambiguous range interval is %g km.rr,rua/1e3);fprin tf(fid,rThere are %g targets with the followi ng paramete

18、rs:, .Ntargets);for i = 1:Ntargetsfprin tf(fid,r ran ge=%5.2g km, SNR=%7.3g dB, rel_RCS=%7.3g dB, vel=%9.4g m/s, .ran ges(i)/1e3,SNR(i),rel_RCS(i),vels(i);endfelose(fid);fprintf(ninData is in file ,data_file)fprin tf(nStude nt data is in file ,mystery_file)fprin tf(nListi ngisinfile,listi ng_file,nn

19、) 到的函数 funetion y = radar( x, fs, T_0, g, T_out, T_ref, fe, r, snr, v)% RADAR simulate radar returns from a si ngle pulse or burst%of ide ntical pulses% usage:% R = radar( X, Fs, T_0, G, T_out, T_ref, Fe, R, SNR, V )% X:baseba nd si ngle pulse waveform (eomplex veetor)%Fs:sampli ng freque ncy of in

20、put pulse in Hz%T_0:start time(s) of in put pulse(s)see%(nu mber of pulses in burst assumed = len gth(g)% G: complex gain(s) of pulse(s),即慢时间，各个 PRI对应的脉冲的前的 加权 20-by-1% T_out: 2-vector T_mi n, T_max defi nes output% T_ref: system refere nee time, n eeded to simulate%burst returns. THIS IS THE t=0 TI

21、ME !% Fc:een ter freq. of the radar.in Hz% R: vector of ran ges to target(s) meters%(nu mber of targets assumed = len gth(r)% SNR: vector of target SNRs (un it no ise power assumed)%This will be SNR *after* allowi ng for RA4% V:vector of target velocities (opti on al)in m/sec%(positive velocities ar

22、e towards the radar)% n ote(1): VELOCITY in meters/sec !%distancesin m, times in sec, BW in Hz.% note(2): assumeseach pulse iscon sta nt (complex) amplitude% no te(3): will accomodate up to quadratic phase pulses% n ote(4): vector of ra nges, R, allows DISTRIBUTED targets%(c) jMcClella n 7/28/90% Mo

23、dified by M. A. Richards, August 1991J = sqrt(-1);c = 3e8;% velocity of light in m/secMx = len gth(x);delta_t = 1/fs;% sampli ng in terval (sec)t_y = T_out(1):delta_t:T_out(2) ;% output sampling times (sec),接收窗的宽度内的等间隔采样337-by-1T_p = Mx*delta_t;% len gth of in put pulse (sec)，基带信号 chirp 的脉冲持续时间，即 Te

24、% Assume zero velocities (stationary targets) if no velocity % vector providedif n argi n 7v = zeros(r);end% en sure that all vectors are colu mn 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 sam

25、plest_x = delta_t*0:(Mx-1);x_ph = unwrap(angle(x); %基带chirp信号的相位，可以看出 x_ph是个抛物线 q = polyfit(t_x,x_ph,2); % 目 的是用 q = axA2 + bx + c 逼近 x_ph% check result using correlation coefficientxfit = polyval(q,t_x); % 看看用 q = axA2 + bx + c 拟合的 xfit 与x_ph的一致程度if (x_ph*xfit)/norm(x_ph)/norm(xfit) = T_out (2) | t

26、max = 0 & t_vals T_p ); %T_p 是 chirp 基带信号的长度，一个 chirp 脉 冲携带有效测量数据,2卩Te上的采样点if tau T_out(2)fpri ntf(nEcho from target #%g at range %g km,i,ri)fprin tf(nFINISHES AFTER the range win dow)non pulse #%g.n,j)end % Place scaled, ra nge-delayed, Doppler shifted pulse into output matrix% Unit no ise power and unit nominal pulse amplitude assumed to% get amplitude from SNR.amp = 10A(s nr(i)/20);% n_out 是对应chirp脉冲宽度Te的120-by-1向量，原来接收的chirp信号未经过脉冲压缩在距 离上占据c*Te/2米的长度，和对应长度的rect信号在距离上占据的长度是一样 的！经过脉冲压缩后chirp信号在距离上才占据c/(2B)米的长度。y(n_out,j) = y(n_out,j) + .( amp * g(j) *exp( -J*2