脉冲压缩,相干积累,MTI,CFAR详解

1、 MATLAB击达数据处理：脉冲压缩，相干积累,M TI,CFAR 详解 2021-05-14 21:42:29| 分类：雷达信号处理|标签：雷达数据仿真 动目标检测 脉压 恒虚警率|字号 订 阅 % Thanks to M. A. Richards for providing the code. clear, hold off format compact J = sqrt(-1); close all % Get root file name for reading results file=input(Enter root file name for data file: ,s); ev

2、al(load ,file,.mat) fprintf(nPulse length = %g microsecondsn,T/1e-6) fprintf(Chirp bandwidth = %g Mhzn,W/1e6) fprintf(Sampling rate = %g Msamples/secn,fs/1e6) figure plot(1e6/fs)*(0:length(s)-1),real(s) imag(s) title(Real and Imaginary Parts of Chirp Pulse) xlabel(time (usec) ylabel(amplitude) grid

3、PRI = 1/PRF; fprintf(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 unam

4、biguous range interval in meters vua = 3e8*PRF/(2*fc); rmin = 3e8*T_out(1)/2; rmax = 3e8*T_out(2)/2; rwin = rmax-rmin; 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) f

5、printf(nThe sampled range window is %g km long.,rwin/1e3) fprintf(nThe unambiguous range interval is %g km.nn,rua/1e3) % Convert range samples to absolute range units. My,Ny=size(y); range=(3e8/2)*(0:My-1)*(1/fs) + T_out(1)/1e3; % T_out(1) 是接收窗的起始时刻，即 range从多远的距离开始观测 pulse = (1:Ny); % Force oversize

6、 FFT, and compute doppler scale factor Lfft = 2A(nextpow2(Ny)+3); % Lfft 是慢时间上 FFT 的点数 doppler = (0:Lfft-1)/Lfft)-0.5)*vua; % 此处的 doppler 表征的是速度，不是 fd = 2v / lambda这个定义式 fprintf(nThe Doppler increment is %g Hz.,PRF/Lfft) fprintf(nThe velocity increment is %g m/s.,3e8*PRF/Lfft/2/fc) % Start with a fe

7、w plots to examine the data % plot power of raw data in dB ydB=db(abs(y)/max(max(abs(y),voltage); figure mesh(pulse,range,ydB) title(FAST-TIME/SLOW-TIME PLOT OF RAW DATA) ylabel(range (km) xlabel(pulse number) % Plot overlay of individual range traces disp() disp() disp(.plotting overlay of range tr

8、aces) figure plot(range,db(y,voltage) title(OVERLAY OF RANGE TRACES) xlabel(distance (km) ylabel(amplitude (dB) grid % Noncoherently integrate the range traces and display disp(.plotting integrated range trace) figure plot(range,db(sum(abs(y).A2),power) % 看看功率的非相干积累 title(NONCOHERENTLY INTEGRATED RA

9、NGE TRACE) xlabel(range bin) ylabel(power) grid % Doppler process and square-law detect the whole % unprocessed array and display mesh. % Use Hamming window throughout. disp(. puting raw range-Doppler map) % 下句效果和 Y=fft(conj(y).*(hamming(Ny)*ones(1,My),Lfft,1); 一样。dim = 1 表示 沿着列做fft Y=fft(conj(y).*(

10、hamming(Ny)*ones(1,My),Lfft); %y 是 20-by-337 的，对 y的慢时间维 (按列)做Lfft点FFT,进行相干积累 % Y= fftshift(Y.*conj(Y),1); 的效果与 Y=git_rotate(Y.*conj(Y),Lfft/2); 完全相同， git_rotate是将矩阵的上下循环移位 Y=git_rotate(Y.*conj(Y),Lfft/2); % note we take mag-squared of Y here also 注意 Y是 做完相干积累后才取模平方的 YdB=db(abs(Y),power); figure mesh

11、(doppler,range,YdB) Since I used the DFT(K points) to perform the Doppler analysis, bin #0 (#1 in MATLAB indexing) corresponds to zero Doppler shift, while each successive bin represents an increment of PRF/K Hz.I used a modest DFT size of 256 to get adequate Doppler definition, the corresponding in

12、crement in velocity is (lambda*PRF/2/K) m/s per FFT bin. The DFT causes the target energy to concentrate along the frequency axis; a target signal that was a sinusoid across all 20 slow time samples and was lost among the noise and clutter is now an asinc function peaking at the appropriate Doppler

13、shift. The asinc is still fairly fat simply because I have only 20 slow-time samples, which limits the Doppler resolution. The noise remains spread out across all of the range and velocity bins. The clutter is now a clearly evident concentration of energy around zero velocity, extending through all

14、of the range bins. RANtE-DOPPLER PLOT OF UNPROCESSED DATA Note also that the four targets are now clearly visible. Each is in the form of a of energy extending across 1.5 km of range. This is because I have not yet done pulse compression. The signal from a target on a given pulse is just an echo of

15、my transmitted pulse, which was 10 肉 long. This is equivalent to 02)(10 s) = 1.5 km. title(RANGE-DOPPLER PLOT OF UNPROCESSED DATA) ylabel(range (km) xlabel(velocity (m/s) levels=(max(YdB(:)+-1 -5 -10 -15 -20 -25 -30); figure contour(doppler,range,YdB,levels) title(RANGE-DOPPLER CONTOUR PLOT OF UNPRO

16、CESSED DATA) ylabel(range (km) xlabel(velocity (m/s) grid % Now start processing the data . % Pulse compression first. Use time-domain Hamming weighting of the % impulse response for range sidelobe control Ls = length(s); % Ls = 120 disp(.performing matched filtering) h = conj( s(Ls:-1:1) ); % 接收信号的

17、时域匹配滤波器：h(t) = conj( s(-t)，注意我以前 不会序列的反转如何实现 h = h.*hamming(length(h); % time-domain Hamming window for range sidelobe control yp = zeros(My+length(h)-1,Ny); % My+length(h)-1 是一个脉冲宽度内的 My 点的 chirp回波信号与h(n)卷积结果的数据长度 for i=1:Ny % Ny是一个CPI中包含的脉冲个数，对每个脉冲做匹配滤波 yp(:,i) = conv(h,y(:,i); end Myp,Nyp=size(yp

18、); % yp 大小为 456-by-20 % yp = fftfilt( h, y ); % using fftfilt instead of conv because it filters %注意yp=fftfilt(h,y) 的结果与yp = conv(h,y)的前length(yp)点相同！此处匹配滤波 要的是conv(h,y)的后面的值 %因此fftfilt在这里并不能得到需要的结果， fftfilt和conv的比拟见最后！ %y1 = fftfilt(b,x); % using fftfilt hump %y2 = conv(b,x,full); y2=y2(1:length(y1

19、); % % multiple columns with one call % compute new range and time scales here to take account of increased range % length due to convolution and offset due to filter delay of Ls-1 samples . %我以前不知道经过长度为 Ls的匹配滤波器 后，输出结果有Ls-1个点的延迟，难怪 计算结果时间上对不上。用固定点数的频域 FFT做匹配滤波不存在延迟这种情况。 不 同于普通长度为 s的FIR滤波器delay = (s

20、-1)/2 sample points ,或者 delay = (s-1)*ts/2 second 。 where,ts refer to sample period and s is tap number. the unit of delay is second. rangep = (3e8/2)*(0:Myp-1)- (Ls-1)*(1/fs) + T_out(1)/1e3; timep = (1e3*rangep)*2/3e8; ypdB=db(abs(yp),voltage); figure mesh(pulse,rangep,ypdB) FAST-TIME/SLOW-TlMf PLO

21、T OF PULSE-COPf?E3SE0 DATA title(FAST-TIME/SLOW-TIME PLOT OF PULSE-COMPRESSED DATA) ylabel(range (km) xlabel(pulse number) levels=(max(ypdB(:)+-1 -5 -10 -15 -20); figure contour(pulse,rangep,ypdB,levels) title(FAST-TIME/SLOW-TIME CONTOUR PLOT OF PULSE-COMPRESSED DATA) ylabel(range (km) xlabel(pulse

22、number) grid % Range-Doppler plots of pulse-compressed data YP=fft(co nj(yp).*(hamming(Nyp)*ones(1,Myp),Lfft); %yp 是 20-b?-456 矩阵,20 代表慢 时间，对该维做Lfft = 256点FFT进行相干积累 YP=git_rotate(YP.*conj(YP),Lfft/2); YPdB=db(abs(YP),power); figure RAMGE-DOPPLER PLOT OF PULSE-COMPRESSED DATA Note that the 1.5 km ext

23、ended humps are now compressed into a peak in a single range bin. The absolute level of the noise, clutter, and target signals are all increased by the pulse compression filtering. The impulse response hasthe same length as the pulse waveform. Linear filtering is a form of integration: the output is

24、 a weighted sum of input samples, with the weights being the impulse response coefficients. For white random processes such as noise and clutter, integrating (summing)samples increases the power by a factor oL. The coherent addition of the target samples by the matched filter increases the target po

25、wer at its peak by2. mesh(doppler,rangep,YPdB) title(RANGE-DOPPLER PLOT OF PULSE-COMPRESSED DATA) ylabel(range (km) xlabel(velocity (m/s) levels=(max(YPdB(:)+-1 -5 -10 -15 -20 -25 -30); figure contour(doppler,rangep,YPdB,levels) title(RANGE-DOPPLER CONTOUR PLOT OF PULSE-COMPRESSED DATA) ylabel(range

26、 (km) xlabel(velocity (m/s) grid % Apply three-pulse canceller in each range bin to raw data disp(.performing 3-pulse clutter cancellation) %yp 是做了脉冲压缩的信号 456-by-20 h = 1 -2 1; % 3 脉冲 MTI 滤波器的 H(z) = 1 - 2z.A(-1) + z.A2 ypm = zeros(Myp,Nyp+length(h)-1); % MTI 滤波器是抑制杂波的，看看 MTI是如何作 用在慢时间上的，即脉冲与脉冲间的， 脉

27、冲间的线性组合原来也要用 conv来完成，虽 然是简单的相加减！ for i=1:Myp ypm(i,:) = conv(h,yp(i,:); % 对于相同的距离单元，在脉冲间做运算：Vout = V(i) - 2V(i -1) + V(i - 2) 是这样干的 end %这里没有去掉 长度为s的FIR滤波器delay = (s-1)/2 sample points(FIR 滤波 器的群时延)对 doppler 频率计算的影响,I did not remove any of the transients at the beginning or end of the convolution, s

28、o the result now has 20+3-1 = 22 samples in slow time. 注意此处3脉冲MTI滤波的结果严格将有误，仅仅去掉滤波结果的前后各 (s-1)/2个采 样点都不够。必须在慢时间维的总脉冲个数 20中减去(s-1=2)才是正确结果！ 20个脉 冲经过3脉冲对消后的结果应该为 20-3+1=18个脉冲。3脉冲对消实际上是hn=1 -2 1与慢时间上的各个脉冲进行脉冲间滤波，即做内积。想想一下，编号为 #1 , #2 , #3 的3个脉冲与hn做内积后变成了一个脉冲，编号为 #2 , #3, #4的3个脉冲与hn 做内积后变成了下一个脉冲.同时应注意hn

29、这个滑动窗必须全部在 20个脉冲之内 才有效，最后#18,#19,#20的脉冲经过3脉冲MTI对消后就结束了。 也就是hn的第 一个系数与1#脉冲重合开始算，到hn的最后一个系数与最后一个脉冲 20#脉冲重合 时结束！对于上述ypm(i,:) = conv(h,yp(i,:); 的结果是22个脉冲，要去掉 ypm(i,:)前 端的(s-1=2)个hn没有完全没入的卷积结果，以及ypm(i,:)后面(s-1=2)个hn没有完 全没入的卷积结果，才是最终正确的剩余 18个脉冲的处理结果。 Mypm,Nypm=size(ypm); % Doppler process and square-law d

30、etect the whole % clutter-cancelled array and display mesh disp(. puting clutter-cancelled range-Doppler map) YPM=fft(conj(ypm).*(hamming(Nypm)*ones(1,Mypm),Lfft); % ypm 是 456-by-22 的， 对慢时间做FFT进行相干积累 YPM=git_rotate(YPM.*conj(YPM),Lfft/2); % 等效 Y= fftshift(YPM.*conj(YPM),1); I have seen many cases wh

31、ere a student does a good job of estimating freque except that they forget to do arfiftshift on their data, or forget that they did do an fftshift on their data. YPMdB=db(abs(YPM),power); figurency, RANGE-DOPPLFt PLOT OF CLUTTEF?-CANCELLED DATA 10DW mesh(doppler,rangep,YPMdB) title(RANGE-DOPPLER PLO

32、T OF CLUTTER-CANCELLED DATA) ylabel(range (km) xlabel(velocity (m/s) levels=(max(YPMdB(:)+-1 -5 -10 -15 -20 -25); figure contour(doppler,rangep,YPMdB,levels) title(RANGE-DOPPLER CONTOUR PLOT OF CLUTTER-CANCELLED DATA) ylabel(range (km) xlabel(velocity (m/s) grid % OK, now lets do a 1D range-only CFA

33、R to set a threshold map. % Define CFAR window 23 in range, including test cell(cell under test) and 2 guard cells(altogether 4 guard cells), so % 18 actual averaging cells cfar = ones(23,1)/18; % 我以前不知道 CFAR是如何实现的 cfar(10:14)=0; % Cell Under Test 单元的 index 是 12 ,它左、右各有 2 个 guard cells , 9 个 averagi

34、ng cells % Now convolve it in range dimension with range-Doppler map YPM % to get average noise map; discard half-transients to maintain data array size N = zeros(size(YPM); for i=1:Lfft temp = conv(YPM(i,:),cfar); %cfar 是 23 点的 h(n) , temp 是长度为 478 的行向量 456 + 23 -1 = 478 N(i,:) = temp(12:end-11); %

35、 用cfar卷积：首先输出会产生 (23-1)/2=11的延迟，即 half-transients ，将其去掉即从temp(12)开始；其次，输出的最后(23-1)/2=11个点是 局部cfar与信号卷积的结果，也是 half-transients ，应该去掉。最终剩下的temp结果 是完整的cfar在信号中的结果，其长度与信号等长， N对噪声功率的估计结果，是 256-by-456的与数据阵 YPM大小相同。 end % plot the resulting noise map NdB=db(abs(N),power); figure mesh(doppler,rangep,NdB) RAN

36、GE-DOPPLER PLOT OF CFAR NOISE ESTIMATE 从上图的结果中可以看出噪声大约为 20dB title(RANGE-DOPPLER PLOT OF CFAR NOISE ESTIMATE) ylabel(range (km) xlabel(velocity (m/s) % multiply the noise by a factor of 16 (20 dB) and then threshold the data *声由率的估计停1是 (7.10) 删要求的门眼W以曲佑计刊的干扰功率栾以 个赃数对到.即 (7.11) 因为T祝凹率-戊相用的翰弱门眼也购冷整单元冏

37、函的临近甲元T扰功率的平均短口的,驴 以以这 BAR方肽称为为尊平均CFARC& CFAR).曲J千拽免率木是精硝巳知而垫通 过估计特到的-旃以k式中的尺度因子口不铮F式02)中的尺度因子，抹将在7J节中避 行推廿 threshold = 16*N; % 我以前不知道20dB是从以上估计出来的 detected = NaN*ones(size(YPM); detected(YPMthreshold) = YPM(YPMthreshold); % plot the resulting detection map on a dB scale figure mesh(doppler,rang

38、ep,db(detected,power) title(RANGE-DOPPLER DETECTION MAP) ylabel(range (km) xlabel(velocity (m/s) % To search for range bins with targets , first noncoherently integrate % across the frequency bins YPMrange = sum(YPM); % 注意YPM是做完相干积累后，取模平方后的结果！， YPMrange是各个距离单元的噪声功率 % noise floor estimate based on me

39、dian to avoid elevation of estimate by % target responses Nrange = median(YPMrange); % YPMrange 是各个脉冲经过非相干积累的结果, Nrange是各个距离单元的噪声功率的中值 Trange = 8*Nrange; % threshold 8x (9 DB) above noise estimate % This loop identifies which range bins have local peaks above the % threshold. It also sets up a vecto

40、r for plotting convenience to circle % the peaks that are found. spikesr = ; % 这一步是寻找距离方向上的峰值 marker_r = NaN*ones(1,length(YPMrange); % 1-by-456 for i=2:Mypm-1 % 为什么是从 YPMrange(2:Mypm-1)? 因为 YPMrange(i)是分别与它 左右两边的值比拟 if ( (YPMrange(i) YPMrange(i+1) & (YPMrange(i) YPMrange(i-1) ) & (YPMrange(

41、i) Trange) spikes r=spikesr; i, YPMrange(i) ; %看看如何给spikesr后面按行续值，这里最 终 spikesr 是 3-by-2 的 marker_r(i) = YPMrange(i); % marker_r 的大局部值是 NaN ,其中有值的是搜索 出的峰值的大小 end end Mspikesr,Nspikesr=size(spikesr); figure plot(rangep,db(YPMrange;Nrange*ones(1,length(YPMrange); . %YPMrange 是非 相干积累结果，Nrange是中值 Trange

42、*ones(1,length(YPMrange),power); %Trange 是阈值 hold on plot(rangep,db(marker_r,power),-ro) % marker_r 中有值的是搜索出的峰值， marker_r 是和YPMrange相同长度的行向量 hold off xlabel(range (km) ylabel(power (dB) title(RANGE PEAKS) grid targets =; % Now find the Doppler peak(s) for each range bin having a target(s). Keep % ad

43、joining Doppler values as well to support subsequent interpolation. for i = 1:Mspikesr % Mspikesr是距离方向搜出的峰值的个数，是 spikes r的行数，这 里 Mspikesr = 3 rb = spikes r(i,1) % current range bin 注意 spikesr 是 3-by-2 的，spikesr 的第一歹U 存的是峰值在距离单元的序号 % search in Doppler in this range bin. Have to do this modulo the % D

44、FT size to allow for Doppler peaks near the spectrum edges. spikes d =; marker_d = NaN*ones(1,Lfft); % 峰值的距离单元上有 Lfft=256个多普勒单元，故 marker_d 是 1-by-256 的 for k=1:Lfft km1 = k - 1; % km1是k左边的那个值 if km1 Lfft kp1 = kp1-Lfft; % modulo the DFT size to allow for Doppler peaks near the spectrum edges end if

45、( (YPM(k,rb) YPM(kp1,rb) & (YPM(k,rb) YPM(km1,rb) ) & (YPM(k,rb) Nrange/2) %Nrange 是 YPMrange 的中值 spikesd=spikesd; k,YPM(km1,rb),YPM(k,rb),YPM(kp1,rb); % k 是峰值所在多 普勒单元的序号，后面 3个值用来2次插值 marker_d(k) = YPM(k,rb); % marker_d 是 1-by-256 的，非 NaN 值对应着多普 勒单元上的峰值大小 end end % end of loop over Doppler b

46、ins Mspikesd,Nspikesd=size(spikesd); figure plot(doppler,db(YPM(:,rb), Nrange*ones(Lfft,1),power); % YPM(:,rb) 是 rb 距离单 元的峰值对应的doppler曲线 hold on plot(doppler,db(marker_d,power),-ro) % marker_d 是搜索出来的 doppler 峰值 hold off xlabel(velocity (m/s) ylabel(power (dB) title(DOPPLER PEAKS FOR RANGE BIN ,int2s

47、tr(rb) grid % For each local peak in Doppler, use peakinterp % to refine the amplitude and location estimate using a quadratic % interpolation for i=1:Mspikesd % note that peakinterp works on magnitude, not magnitude-squared amp del_k = peakinterp(sqrt(spikesd(i,2:4); targets = targets;amp, rangep(r

48、b), spikesd(i,1)+del_k; end end % end of loop over range bins containing peaks Mtargets, Ntarget s = size(targets); fprintf(nnINTERIM PARAMETERS OF DETECTED TARGETS:n) fprintf(nNumber Power Range (km) DFT Index) for i = 1:Mtargets fprintf(n %2.0g %7.3g %6.3g %9.4g, . i,targets(i,:) end disp() disp()

49、 % OK, now targets contains a list, we hope, of all targets with the % range bin, interpolated peak magnitude, and interpolated Doppler bin of % each. Now we need to get in the desired units and adjust the amplitudes % for a couple of factors. % Put Doppler peaks into m/s units DFT indices k are con

50、verted to velocity using the relationship v = k( /2)(PRF/K). Remember thatk counts from 0 in this formula, whereas MATLAB counts from 1. targets(:,3) = (targets(:,3)-1)/Lfft)-0.5)*vua; % Adjust target amplitudes for MTI filter, % and then convert to relative RCS targets(:,1) = targets(:,1)./(4*(sin(

51、targets(:,3)*(pi/vua).A2); fprintf(nnINTERIM PARAMETERS OF DETECTED TARGETS AFTER MTI FILTER CORRECTION) fprintf(nNumber Power Range (km) Vel (m/s) for i = 1:Mtargets fprintf(n %2.0g %7.3g %6.3g %9.4g, . i,targets(i,:) end disp() disp() Note that theR4 is the correct scaling for target data, but not

52、 for clutter, which would require either R2 or R3, depending on whether it was beam- or pulse-limited clutter. (Earlier I noted that the code that generated the data assumed beam-limited but you don t know that based on just the information provided to you.) However, I have filtered out most of the

53、clutter and I am concerned with target amplitudes, not clutter amplitude. Also, I applied the range correction to power (magnitude-squared) data. If I wanted to do the correction on voltage data, I would usR2 instead of R4. targets(:,1) = targets(:,1).*targets(:,2).A2; targets(:,1) = targets(:,1)/ma

54、x(targets(:,1); targets(:,1) = db(targets(:,1),voltage); % List out detected target amplitudes, ranges, Dopplers fprintf(nnESTIMATED PARAMETERS OF DETECTED TARGETS:n) fprintf(nNumber Rel RCS (dB) Range (km) Vel (m/s) for i = 1:Mtargets fprintf(n %2.0g %7.3g %6.3g %9.4g, . i,targets(i,:) end disp() disp() 用到的函数 function amp,del_index = peakinterp(z) %PEAKINTERP % peakinterp performs a quadratic interpolation of the % peak