Matlab雷达回波数据模拟

1、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;      % ch

2、irp 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); % 

3、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);  &

4、#160;    % pulse index array, 慢时间采样的序列，注意第一个PRI标记为0是为了慢时间起始时刻从零开始PRF = 10.0e3;         % PRF in HzPRI = (1/PRF);        % PRI in secT_0 = PRI*jkl;        % relative start time

5、s 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_a

6、t_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.

7、9;,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

8、= 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,

9、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;                

10、60;  % 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;  

11、        % 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 ar

12、e %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 b

13、in - 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

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

15、;. 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) );% For

16、ce 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 =

17、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

18、 = 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

19、 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(&#

20、39;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 ba

21、ndwidth = %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 limi

22、ts 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,'

23、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', .     

24、;      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

25、, 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 )%

26、      X:      baseband single pulse waveform (complex vector)%      Fs:     sampling frequency of input pulse        in Hz%      T_0:  

27、  start time(s) of input pulse(s)          sec%                (number of pulses in burst assumed = length(g) )%      G:     

28、 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%     

29、; T_ref:  system "reference" time, needed to simulate%                 burst returns. THIS IS THE "t=0" TIME !%      Fc:     center freq. of the radar

30、.                in Hz%      R:      vector of ranges to target(s)             meters%     

31、;            (number of targets assumed = length(r) )%    SNR:      vector of target SNRs (unit noise power assumed)%              This will

32、 be SNR *after* allowing for R4%      V:      vector of target velocities (optional)    in m/sec%                 (positive velocities are towards the rada

33、r)%  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 r

34、anges, R, allows DISTRIBUTED targets%     (c) jMcClellan 7/28/90%         Modified by M. A. Richards, August 1991J = sqrt(-1);c = 3e8;                 

35、0;             % velocity of light in m/secMx = length(x);delta_t = 1/fs;                        % sampling interval (sec)t_y

36、 = 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

37、 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% lat

38、er 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_p

39、h)/norm(xfit) < 0.99    disp('pulse phase is not quadratic')    keyboardend%-  Form (initially empty) output matrix -%Mr = length(t_y);  Nr = length(g);     % output samples in a matrixy = zeros(Mr,Nr); % 337-by-20% Index 'i&

40、#39; loops over the number of targetsfor  i = 1:length(r)    ri = r(i);    vi = v(i);    f_doppler = 2*vi*fc/c;  %看看目标以速度v接近雷达应该如何表示，f_doppler = 2*v/lambda    % Index 'j' loops over the number of pulses    f

41、or j = 1:length(g) % 20 pulses        r_at_T_0 = ri - vi*T_0(j); %一般 T_0 = 0表示当下目标离雷达距ri，若 T_0 = 7 则表示7秒前目标离雷达距ri，现在目标距雷达？        % Compute start and end time of reflected pulse at receiver,     

42、   % ensure that it falls at least partially within the range (time) window        tau = 2*r_at_T_0/(c+vi);   tmax = tau + T_p;        if tau >= T_out(2) | tmax <= T_out(1)     &#

43、160;      fprintf('nEcho from target #%g at range %g km',i,ri)            fprintf('nis COMPLETELY OUT OF the range window')            fprintf('n

44、on pulse #%g.n',j)        else            % Figure out which sample locations in the output grid contain            % reflected pulse   

45、         t_vals = t_y - tau;  %t_y是接收窗的区间  tau = 2*r_at_T_0/(c+vi) ，表示tau时刻目标位于r_at_T_0 ，tau即书上的t0 = 2R0/c            % 用t_vals 而不是t_y是为了用 t_vals与 0, T_p这两个上下界相比较，可以想象若用t_y应该和 tau 和 T_p

46、 + tau比较            n_out = find( t_vals >= 0  &  t_vals < T_p ); %T_p是chirp基带信号的长度，一个chirp脉冲携带有效测量数据，即Te上的采样点            if tau < T_out(1)  &

47、#160;             fprintf('nEcho from target #%g at range %g km',i,ri)                fprintf('nSTARTS BEFORE the range window')   

48、60;            fprintf('non pulse #%g.n',j)            end            if tmax > T_out(2)     &

49、#160;          fprintf('nEcho from target #%g at range %g km',i,ri)                fprintf('nFINISHES AFTER the range window')      &#

50、160;         fprintf('non pulse #%g.n',j)            end            % Place scaled, range-delayed, Doppler shifted pulse into output matr

51、ix            % Unit noise power and unit nominal pulse amplitude assumed to            % get amplitude from SNR.            amp = 10(s

52、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) + .                                 &#