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

MATLAB 雷达数据处理 脉冲压缩 相干积累 M TI CFAR 详解 2014 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 eval load file mat fprintf 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 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 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 unambiguous 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 fprintf nThe sampled range window is g km long rwin 1e3 fprintf nThe unambiguous range interval is g km n n 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 FFT and compute doppler scale factor Lfft 2 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 few 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 traces 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 2 power 看看功率的非相干积累 title NONCOHERENTLY INTEGRATED RANGE 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 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 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 increment 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 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 of the range bins Note also that the four targets are now clearly visible Each is in the form of a hump 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 my transmitted pulse which was 10 s long This is equivalent to c 2 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 UNPROCESSED 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 接收信号的时域匹配滤波器 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 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 y2 conv b x full y2 y2 1 length y1 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 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 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 number grid Range Doppler plots of pulse compressed data YP fft conj yp hamming Nyp ones 1 Myp Lfft yp 是 20 by 456 矩阵 20 代表 慢时间 对该维做 Lfft 256 点 FFT 进行相干积累 YP git rotate YP conj YP Lfft 2 YPdB db abs YP power figure Note that the 1 5 km extended 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 has the same length as the pulse waveform Linear filtering is a form of integration the output is 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 L samples increases the power by a factor of L The coherent addition of the target samples by the matched filter increases the target power at its peak by L2 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 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 1 z 2 ypm zeros Myp Nyp length h 1 MTI 滤波器是抑制杂波的 看看 MTI 是如何作 用在慢时间上的 即脉冲与脉冲间的 脉冲间的线性组合原来也要用脉冲间的线性组合原来也要用 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 so 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 脉冲对消实际上是脉冲对消实际上是 h n 1 2 1 与慢时间上的各个脉冲进行脉冲间滤波 即做内积 想想一下 编号为与慢时间上的各个脉冲进行脉冲间滤波 即做内积 想想一下 编号为 1 2 3 的的 3 个脉冲与个脉冲与 h n 做内积后变成了一个脉冲 编号为做内积后变成了一个脉冲 编号为 2 3 4 的的 3 个个 脉冲与脉冲与 h n 做内积后变成了下一个脉冲做内积后变成了下一个脉冲 同时应注意同时应注意 h n 这个滑动窗必须全部在这个滑动窗必须全部在 20 个脉冲之内才有效 最后个脉冲之内才有效 最后 18 19 20 的脉冲经过的脉冲经过 3 脉冲脉冲 MTI 对消后就结束了 对消后就结束了 也就也就 是是 h n 的第一个系数与的第一个系数与 1 脉冲重合开始算 到脉冲重合开始算 到 h n 的最后一个系数与最后一个脉冲的最后一个系数与最后一个脉冲 20 脉冲重合时结束 脉冲重合时结束 对于上述对于上述 ypm i conv h yp i 的结果是的结果是 22 个脉冲 要去个脉冲 要去 掉掉 ypm i 前端的前端的 s 1 2 个个 h n 没有完全没入的卷积结果 以及没有完全没入的卷积结果 以及 ypm i 后面后面 s 1 2 个个 h n 没有完全没入的卷积结果 才是最终正确的剩余没有完全没入的卷积结果 才是最终正确的剩余 18 个脉冲的处理结果 个脉冲的处理结果 Mypm Nypm size ypm Doppler process and square law detect 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 where a student does a good job of estimating frequency except that they forget to do an fftshift on their data or forget that they did do an fftshift on their data YPMdB db abs YPM power figure mesh doppler rangep YPMdB title RANGE DOPPLER PLOT 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 let s do a 1D range only CFAR 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 个 averaging 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 用 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 从上图的结果中可以看出噪声大约为从上图的结果中可以看出噪声大约为 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 threshold 16 N 我以前不知道 20dB 是从以上估计出来的 detected NaN ones size YPM detected YPM threshold YPM YPM threshold plot the resulting detection map on a dB scale figure mesh doppler rangep 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 median 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 vector 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 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 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 adjoining Doppler values as well to support subsequent interpolation for i 1 Mspikesr Mspikesr 是距离方向搜出的峰值的个数 是 spikesr 的行数 这 里 Mspikesr 3 rb spikesr i 1 current range bin 注意 spikesr 是 3 by 2 的 spikesr 的第一列 存的是峰值在距离单元的序号 search in Doppler in this range bin Have to do this modulo the DFT size to allow for Doppler peaks near the spectrum edges spikesd 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 YPM k rb YPM kp1 rb 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 bins 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 int2str 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 rb spikesd i 1 del k end end end of loop over range bins containing peaks Mtargets Ntargets size targets fprintf n nINTERIM 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 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 converted to velocity using the relationship v k 2 PRF K Remember that k 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 targets 3 pi vua 2 fprintf n nINTERIM PARAMETERS OF DETECTED TARGETS AFTER MTI FILTER CORRECTION n 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 the R4 is the correct scaling for target data but not 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 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 useR2 instead of R4 targets 1 targets 1 targets 2 2 targets 1 targets 1 max targets 1 targets 1 db targets 1 voltage List out detected target amplitudes ranges Dopplers fprintf n nESTIMATED PARAMETERS OF DETECTED TARGETS n fprintf nNumber Rel RCS dB Range km Vel m s for i 1 Mtargets fprintf n 2 0g 7