脈沖壓縮,相干積累,MTI,CFAR詳解

1、 MATLAB擊達(dá)數(shù)據(jù)處理：脈沖壓縮，相干積累,M TI,CFAR 詳解 2021-05-14 21:42:29| 分類：雷達(dá)信號(hào)處理|標(biāo)簽：雷達(dá)數(shù)據(jù)仿真 動(dòng)目標(biāo)檢測 脈壓 恒虛警率|字號(hào) 訂 閱 % 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) 是接收窗的起始時(shí)刻，即 range從多遠(yuǎn)的距離開始觀測 pulse = (1:Ny); % Force oversize

6、 FFT, and compute doppler scale factor Lfft = 2A(nextpow2(Ny)+3); % Lfft 是慢時(shí)間上 FFT 的點(diǎn)數(shù) doppler = (0:Lfft-1)/Lfft)-0.5)*vua; % 此處的 doppler 表征的是速度，不是 fd = 2v / lambda這個(gè)定義式 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的慢時(shí)間維 (按列)做Lfft點(diǎn)FFT,進(jìn)行相干積累 % Y= fftshift(Y.*conj(Y),1); 的效果與 Y=git_rotate(Y.*conj(Y),Lfft/2); 完全相同， git_rotate是將矩陣的上下循環(huán)移位 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) ); % 接收信號(hào)的

17、時(shí)域匹配濾波器：h(t) = conj( s(-t)，注意我以前 不會(huì)序列的反轉(zhuǎn)如何實(shí)現(xiàn) h = h.*hamming(length(h); % time-domain Hamming window for range sidelobe control yp = zeros(My+length(h)-1,Ny); % My+length(h)-1 是一個(gè)脈沖寬度內(nèi)的 My 點(diǎn)的 chirp回波信號(hào)與h(n)卷積結(jié)果的數(shù)據(jù)長度 for i=1:Ny % Ny是一個(gè)CPI中包含的脈沖個(gè)數(shù)，對每個(gè)脈沖做匹配濾波 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) 的結(jié)果與yp = conv(h,y)的前l(fā)ength(yp)點(diǎn)相同！此處匹配濾波 要的是conv(h,y)的后面的值 %因此fftfilt在這里并不能得到需要的結(jié)果， 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 . %我以前不知道經(jīng)過長度為 Ls的匹配濾波器 后，輸出結(jié)果有Ls-1個(gè)點(diǎn)的延遲，難怪 計(jì)算結(jié)果時(shí)間上對不上。用固定點(diǎn)數(shù)的頻域 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 代表慢 時(shí)間，對該維做Lfft = 256點(diǎn)FFT進(jìn)行相干積累 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 是做了脈沖壓縮的信號(hào) 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是如何作 用在慢時(shí)間上的，即脈沖與脈沖間的， 脈

27、沖間的線性組合原來也要用 conv來完成，雖 然是簡單的相加減！ for i=1:Myp ypm(i,:) = conv(h,yp(i,:); % 對于相同的距離單元，在脈沖間做運(yùn)算：Vout = V(i) - 2V(i -1) + V(i - 2) 是這樣干的 end %這里沒有去掉 長度為s的FIR濾波器delay = (s-1)/2 sample points(FIR 濾波 器的群時(shí)延)對 doppler 頻率計(jì)算的影響,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濾波的結(jié)果嚴(yán)格將有誤，僅僅去掉濾波結(jié)果的前后各 (s-1)/2個(gè)采 樣點(diǎn)都不夠。必須在慢時(shí)間維的總脈沖個(gè)數(shù) 20中減去(s-1=2)才是正確結(jié)果！ 20個(gè)脈 沖經(jīng)過3脈沖對消后的結(jié)果應(yīng)該為 20-3+1=18個(gè)脈沖。3脈沖對消實(shí)際上是hn=1 -2 1與慢時(shí)間上的各個(gè)脈沖進(jìn)行脈沖間濾波，即做內(nèi)積。想想一下，編號(hào)為 #1 , #2 , #3 的3個(gè)脈沖與hn做內(nèi)積后變成了一個(gè)脈沖，編號(hào)為 #2 , #3, #4的3個(gè)脈沖與hn 做內(nèi)積后變成了下一個(gè)脈沖.同時(shí)應(yīng)注意hn

29、這個(gè)滑動(dòng)窗必須全部在 20個(gè)脈沖之內(nèi) 才有效，最后#18,#19,#20的脈沖經(jīng)過3脈沖MTI對消后就結(jié)束了。 也就是hn的第 一個(gè)系數(shù)與1#脈沖重合開始算，到hn的最后一個(gè)系數(shù)與最后一個(gè)脈沖 20#脈沖重合 時(shí)結(jié)束！對于上述ypm(i,:) = conv(h,yp(i,:); 的結(jié)果是22個(gè)脈沖，要去掉 ypm(i,:)前 端的(s-1=2)個(gè)hn沒有完全沒入的卷積結(jié)果，以及ypm(i,:)后面(s-1=2)個(gè)hn沒有完 全沒入的卷積結(jié)果，才是最終正確的剩余 18個(gè)脈沖的處理結(jié)果。 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 的， 對慢時(shí)間做FFT進(jìn)行相干積累 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是如何實(shí)現(xiàn)的 cfar(10:14)=0; % Cell Under Test 單元的 index 是 12 ,它左、右各有 2 個(gè) guard cells , 9 個(gè) 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 點(diǎn)的 h(n) , temp 是長度為 478 的行向量 456 + 23 -1 = 478 N(i,:) = temp(12:end-11); %

35、 用cfar卷積：首先輸出會(huì)產(chǎn)生 (23-1)/2=11的延遲，即 half-transients ，將其去掉即從temp(12)開始；其次，輸出的最后(23-1)/2=11個(gè)點(diǎn)是 局部cfar與信號(hào)卷積的結(jié)果，也是 half-transients ，應(yīng)該去掉。最終剩下的temp結(jié)果 是完整的cfar在信號(hào)中的結(jié)果，其長度與信號(hào)等長， N對噪聲功率的估計(jì)結(jié)果，是 256-by-456的與數(shù)據(jù)陣 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 從上圖的結(jié)果中可以看出噪聲大約為 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 *聲由率的估計(jì)停1是 (7.10) 刪要求的門眼W以曲佑計(jì)刊的干擾功率欒以 個(gè)贓數(shù)對到.即 (7.11) 因?yàn)門祝凹率-戊相用的翰弱門眼也購冷整單元冏

37、函的臨近甲元T擾功率的平均短口的,驢 以以這 BAR方肽稱為為尊平均CFARC& CFAR).曲J千拽免率木是精硝巳知而墊通 過估計(jì)特到的-旃以k式中的尺度因子口不錚F式02)中的尺度因子，抹將在7J節(jié)中避 行推廿 threshold = 16*N; % 我以前不知道20dB是從以上估計(jì)出來的 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是做完相干積累后，取模平方后的結(jié)果！， YPMrange是各個(gè)距離單元的噪聲功率 % noise floor estimate based on me

39、dian to avoid elevation of estimate by % target responses Nrange = median(YPMrange); % YPMrange 是各個(gè)脈沖經(jīng)過非相干積累的結(jié)果, Nrange是各個(gè)距離單元的噪聲功率的中值 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)? 因?yàn)?YPMrange(i)是分別與它 左右兩邊的值比擬 if ( (YPMrange(i) YPMrange(i+1) & (YPMrange(i) YPMrange(i-1) ) & (YPMrange(

41、i) Trange) spikes r=spikesr; i, YPMrange(i) ; %看看如何給spikesr后面按行續(xù)值，這里最 終 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 是非 相干積累結(jié)果，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是距離方向搜出的峰值的個(gè)數(shù)，是 spikes r的行數(shù)，這 里 Mspikesr = 3 rb = spikes r(i,1) % current range bin 注意 spikesr 是 3-by-2 的，spikesr 的第一歹U 存的是峰值在距離單元的序號(hào) % 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個(gè)多普勒單元，故 marker_d 是 1-by-256 的 for k=1:Lfft km1 = k - 1; % km1是k左邊的那個(gè)值 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 是峰值所在多 普勒單元的序號(hào)，后面 3個(gè)值用來2次插值 marker_d(k) = YPM(k,rb); % marker_d 是 1-by-256 的，非 NaN 值對應(yīng)著多普 勒單元上的峰值大小 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 距離單 元的峰值對應(yīng)的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() 用到的函數(shù) function amp,del_index = peakinterp(z) %PEAKINTERP % peakinterp performs a quadratic interpolation of the % peak