Matlab雷达回波数据模拟

Matlab雷达回波数据模拟

clear, hold off format compact J = 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, seconds W = 10e6; % chirp bandwidth, Hz

fs = 12e6; % chirp sampling rate, Hz; oversample by a littlefprintf('\\nPulse length = %g microseconds\\n',T/1e-6)

fprintf('Chirp bandwidth = %g Mhz\\n',W/1e6)

fprintf('Sampling rate = %g Msamples/sec\\n',fs/1e6) s = git_chirp(T,W,fs/W); % 120-by-1 array plot((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 pulses

jkl = 0:(Np-1); % pulse index array, 慢时间采样的序列，注意第一个PRI标记为0是为了慢时间起始时刻从零开始 PRF = 10.0e3; % PRF in Hz PRI = (1/PRF); % PRI in sec

T_0 = PRI*jkl; % relative start times of pulses, in sec g = ones(1,Np); % gains of pulses

T_out = [12 40]*1e-6; % start and end times of range window in sec， 这个就是接收窗的时间宽度Trec

T_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.\\n\\n',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 km

% (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 R^4

% 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/sec Mx = 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-1

T_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 = ax^2 + bx + c 逼近x_ph% check result using correlation coefficientxfit = polyval(q,t_x); % 看看用 q = ax^2 + bx + c 拟合的xfit与x_ph的一致程度

if (x_ph'*xfit)/norm(x_ph)/norm(xfit) < 0.99 disp('pulse phase is not quadratic') keyboard end%

%--- Form (initially empty) output matrix --- %

Mr = length(t_y); Nr = length(g); % output samples in a matrix

y = zeros(Mr,Nr); % 337-by-20% Index 'i' loops over the number of targets for 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 for 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,

% 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)

fprintf('\\nEcho from target #%g at range %g km',i,ri) fprintf('\\nis COMPLETELY OUT OF the range window') fprintf('\\non pulse #%g.\\n',j)

else % Figure out which sample locations in the output grid contain % reflected pulse

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

fprintf('\\nEcho from target #%g at range %g km',i,ri) fprintf('\\nSTARTS BEFORE the range window') fprintf('\\non pulse #%g.\\n',j) end

if tmax > 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,j) = y(n_out,j) + ... ( amp * g(j) *

exp( -J*2*pi*fc*tau ) ) ... % fc*tau是f0*R0，载频f0乘以R0，是下图的exp中的第一项 .* [ exp( J*2*pi*f_doppler*t_y(n_out) ) ] ... % t_y是接收窗的范围，是下图exp中的第二项 .* [ exp( J*polyval(q,t_vals(n_out)) ) ]; %数字化的chirp信号不能用具体表达式来描述回波的相位，如此给出接收窗t_vals内一个脉冲Te对应的相位，恰好是下图中s(t-t0)的那个复包络！

end % end of \"if

tau >= T_out(2) ....\"

end % end of loop over j

end % end of loop over i function x = git_chirp( T, W, p ) %CHIRP generate a sampled chirp signal % X = git_chirp( T, W,

)

% X: N=pTW samples of a \"chirp\" signal % exp(j(W/T)pi*t^2) -T/2 <= t < +T/2 % T: time duration from -T/2 to +T/2

% W: swept bandwidth from -W/2 to +W/2 % optional:

% P: samples at P times the Nyquist rate (W) % i.e., sampling interval is 1/(PW) % default is P = 1 %

if nargin < 3 p = 1; end

J = sqrt(-1); %-------------- delta_t = 1/(p*W);

N = round( p*T*W ); %--same as T/delta_t, but rounded nn = [0:N-1]';

x = exp( J*pi*W/T * (delta_t*nn - (N-1)/2/p/W).^2 ); % symmetric version