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

1、MATLAB雷達(dá)數(shù)據(jù)處理：脈沖壓縮，相干積累,MTI,CFAR詳解2014-05-1421:42:29|分類：雷達(dá)信號處理|標(biāo)簽：雷達(dá)數(shù)據(jù)仿真動目標(biāo)檢測脈壓恒虛警率|字號訂閱%ThankstoM.A.Richardsforprovidingthecode.clear,holdoffformatcompactJ=sqrt(-1);closeall%Getrootfilenameforreadingresultsfile=input('Enterrootfilenamefordatafile:','s');eval('load',file,'

2、.mat')fprintf('nPulselength=%gmicrosecondsn',T/1e-6)fprintf('Chirpbandwidth=%gMhzn',W/1e6)fprintf('Samplingrate=%gMsamples/secn',fs/1e6)figureplot(1e6/fs)*(0:length(s)-1),real(s)imag(s)title('RealandImaginaryPartsofChirpPulse')xlabel('time(usec)')ylabel(&#

3、39;amplitude')gridPRI=1/PRF;fprintf('nWearesimulating%gpulsesatanRFof%gGHz',Np,fc/1e9)fprintf('nandaPRFof%gkHz,givingaPRIof%gusec.',PRF/1e3,PRI/1e-6)fprintf('nTherangewindowlimitsare%gto%gusec.n',.T_out(1)/1e-6,T_out(2)/1e-6)%ComputeunambiguousDopplerintervalinm/sec%Compu

4、teunambiguousrangeintervalinmetersvua=3e8*PRF/(2*fc);rmin=3e8*T_out(1)/2;rmax=3e8*T_out(2)/2;rwin=rmax-rmin;rua=3e8/2/PRF;fprintf('nTheunambiguousvelocityintervalis%gm/s.',vua)fprintf('nTherangewindowstartsat%gkm.',rmin/1e3)fprintf('nTherangewindowendsat%gkm.',rmax/1e3)fprint

5、f('nThesampledrangewindowis%gkmlong.',rwin/1e3)fprintf('nTheunambiguousrangeintervalis%gkm.nn',rua/1e3)%Convertrangesamplestoabsoluterangeunits.My,Ny=size(y);range=(3e8/2)*(0:My-1)*(1/fs)+T_out(1)/1e3;%T_out(1)是接收窗的起始時刻，即range從多遠(yuǎn)的距離開始觀測pulse=(1:Ny);%ForceoversizeFFT,andcomputedoppler

6、scalefactorLfft=2A(nextpow2(Ny)+3);%Lfft是慢時間上FFT的點(diǎn)數(shù)doppler=(0:Lfft-1)/Lfft)-0.5)*vua;%此處的doppler表征的是速度，不是fd=2v/lambda這個定義式fprintf('nTheDopplerincrementis%gHz.',PRF/Lfft)fprintf('nThevelocityincrementis%gm/s.',3e8*PRF/Lfft/2/fc)%Startwithafewplotstoexaminethedata%plotpowerofrawdataind

7、BydB=db(abs(y)/max(max(abs(y),'voltage');figuremesh(pulse,range,ydB)title('FAST-TIME/SLOW-TIMEPLOTOFRAWDATA')ylabel('range(km)')xlabel('pulsenumber')%Plotoverlayofindividualrangetracesdisp('')disp('')disp('plottingoverlayofrangetraces')figurepl

8、ot(range,db(y,'voltage')title('OVERLAYOFRANGETRACES')xlabel('distance(km)')ylabel('amplitude(dB)')grid%Noncoherentlyintegratetherangetracesanddisplaydisp('plottingintegratedrangetrace')figureplot(range,db(sum(abs(y).A2)')','power')%看看功率的非相干積累ti

9、tle('NONCOHERENTLYINTEGRATEDRANGETRACE')xlabel('rangebin')ylabel('power')grid%Dopplerprocessandsquare-lawdetectthewhole%unprocessedarrayanddisplaymesh.%UseHammingwindowthroughout.disp('computingrawrange-Dopplermap')%下句效果和Y=fft(conj(y').*(hamming(Ny)*ones(1,My),Lff

10、t,1);一樣。dim=1表示沿著列做fftY=fft(conj(y').*(hamming(Ny)*ones(1,My),Lfft);%y是20-by-337的，對y'的慢時間維(按列)做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);%notewetakemag-squaredofYherealso注意Y是做完相干積累后才取模平方的YdB=db(abs(Y

11、),'power');figuremesh(doppler,range,YdB')SinceIusedtheDFT(Kpoints)toperformtheDoppleranalysis,bin#0(#1inMATLABindexing)correspondstozeroDopplershift,whileeachsuccessivebinrepresentsanincrementofPRF/KHz.IusedamodestDFTsizeof256togetadequateDopplerdefinition,thecorrespondingincrementinvelo

12、cityis(lambda*PRF/2/K)m/sperFFTbin.TheDFTcausesthetargetenergytoconcentratealongthefrequencyaxis;atargetsignalthatwasasinusoidacrossall20slowtimesamplesandwaslostamongthenoiseandclutterisnowanasincfunctionpeakingattheappropriateDopplershift.Theasincisstillfairly"fat"simplybecauseIhaveonly2

13、0slow-timesamples,whichlimitstheDopplerresolution.Thenoiseremainsspreadoutacrossalloftherangeandvelocitybins.Theclutterisnowaclearlyevidentconcentrationofenergyaroundzerovelocity,extendingthroughalloftherangebins.RANt?E-DOPPLERPLOTOFUNPROCESSEDDATAhump'Notealsothatthefourtargetsarenowclearlyvisi

14、ble.Eachisintheformofa'ofenergyextendingacross1.5kmofrange.ThisisbecauseIhavenotyetdonepulsecompression.Thesignalfromatargetonagivenpulseisjustanechoofmytransmittedpulse,whichwas10%long.Thisisequivalentto(V2)(10s)=1.5km.title('RANGE-DOPPLERPLOTOFUNPROCESSEDDATA')ylabel('range(km)'

15、;)xlabel('velocity(m/s)')levels=(max(YdB(:)+-1-5-10-15-20-25-30);figurecontour(doppler,range,YdB',levels)title('RANGE-DOPPLERCONTOURPLOTOFUNPROCESSEDDATA')ylabel('range(km)')xlabel('velocity(m/s)')grid%Nowstartprocessingthedata.%Pulsecompressionfirst.Usetime-domai

16、nHammingweightingofthe%impulseresponseforrangesidelobecontrolLs=length(s);%Ls=120disp('performingmatchedfiltering')h=conj(s(Ls:-1:1);%接收信號的時域匹配濾波器：h(t)=conj(s(-t)，注意我以前不會序列的反轉(zhuǎn)如何實(shí)現(xiàn)h=h.*hamming(length(h);%time-domainHammingwindowforrangesidelobecontrolyp=zeros(My+length(h)-1,Ny);%My+length(h)-

17、1是一個脈沖寬度內(nèi)的My點(diǎn)的chirp回波信號與h(n)卷積結(jié)果的數(shù)據(jù)長度fori=1:Ny%Ny是一個CPI中包含的脈沖個數(shù)，對每個脈沖做匹配濾波yp(：,i)=conv(h,y(:,i);endMyp,Nyp=size(yp);%yp大小為456-by-20%yp=fftfilt(h,y);%usingfftfiltinsteadofconvbecauseitfilters%注意yp=fftfilt(h,y)的結(jié)果與yp=conv(h,y)的前l(fā)ength(yp)點(diǎn)相同！此處匹配濾波要的是conv(h,y)的后面的值%因此fftfilt在這里并不能得到需要的結(jié)果，fftfilt和conv的

18、比較見最后！%y1=fftfilt(b,x);%usingfftfilt%y2=conv(b,x,'full');y2=y2(1:length(y1);%multiplecolumnswithonecall%computenewrangeandtimescalesheretotakeaccountofincreasedrange%lengthduetoconvolutionandoffsetduetofilterdelayofLs-1samples.%我以前不知道經(jīng)過長度為Ls的匹配濾波器后，輸出結(jié)果有Ls-1個點(diǎn)的延遲，難怪計(jì)算結(jié)果時間上對不上。用固定點(diǎn)數(shù)的頻域FFT做匹配濾

19、波不存在延遲這種情況。不同于普通長度為s的FIR濾波器delay=(s-1)/2samplepoints,或者delay=(s-1)*ts/2second。where,tsrefertosampleperiodandsistapnumber.theunitofdelayissecond.rangep=(3e8/2)*(0:Myp-1)-(Ls-1)*(1/fs)+T_out(1)/1e3;timep=(1e3*rangep)*2/3e8;ypdB=db(abs(yp),'voltage');figuremesh(pulse,rangep,ypdB)FAST-TIME/SLOW-

20、TlMfPLOTOFPULSE-COPf?E3SE0DATAtitle('FAST-TIME/SLOW-TIMEPLOTOFPULSE-COMPRESSEDDATA')ylabel('range(km)')xlabel('pulsenumber')levels=(max(ypdB(:)+-1-5-10-15-20);figurecontour(pulse,rangep,ypdB,levels)title('FAST-TIME/SLOW-TIMECONTOURPLOTOFPULSE-COMPRESSEDDATA')ylabel(&#

21、39;range(km)')xlabel('pulsenumber')grid%Range-Dopplerplotsofpulse-compresseddataYP=fft(conj(yp').*(hamming(Nyp)*ones(1,Myp),Lfft);%yp是20-by-456矩陣,20代表慢時間，對該維做Lfft=256點(diǎn)FFT進(jìn)行相干積累YP=git_rotate(YP.*conj(YP),Lfft/2);YPdB=db(abs(YP),'power');figureRAMGE-DOPPLERPLOTOFPULSE-COMPRESSE

22、DDATAWONotethatthe1.5kmextendedhumpsarenowcompressedintoapeakinasinglerangebin.Theabsolutelevelofthenoise,clutter,andtargetsignalsareallincreasedbythepulsecompressionfiltering.Theimpulseresponsehasthesamelengthasthepulsewaveform.Linearfilteringisaformofintegration:theoutputisaweightedsumofinputsampl

23、es,withtheweightsbeingtheimpulseresponsecoefficients.Forwhiterandomprocessessuchasnoiseandclutter,integrating(summing)samplesincreasesthepowerbyafactoroL.Thecoherentadditionofthetargetsamplesbythematchedfilterincreasesthetargetpoweratitspeakby2.mesh(doppler,rangep,YPdB')title('RANGE-DOPPLERP

24、LOTOFPULSE-COMPRESSEDDATA')ylabel('range(km)')xlabel('velocity(m/s)')levels=(max(YPdB(:)+-1-5-10-15-20-25-30);figurecontour(doppler,rangep,YPdB',levels)title('RANGE-DOPPLERCONTOURPLOTOFPULSE-COMPRESSEDDATA')ylabel('range(km)')xlabel('velocity(m/s)')gri

25、d%Applythree-pulsecancellerineachrangebintorawdatadisp('.performing3-pulsecluttercancellation')%yp是做了脈沖壓縮的信號456-by-20h=1-21'%3脈沖MTI濾波器的H(z)=1-2z.A(-1)+z.=ypm=zeros(Myp,Nyp+length(h)-1);%MTI濾波器是抑制雜波的，看看MTI是如何作用在慢時間上的，即脈沖與脈沖間的，脈沖間的線性組合原來也要用conv來完成，雖然是簡單的相加減！fori=1:Mypypm(i,:)=conv(h,yp(i,:

26、);%對于相同的距離單元，在脈沖間彳&運(yùn)算：Vout=V(i)-2V(i-1)+V(i-2)是這樣干的end%這里沒有去掉長度為s的FIR濾波器delay=(s-1)/2samplepoints(FIR濾波器的群時延)對doppler頻率計(jì)算的影響,Ididnotremoveanyofthetransientsatthebeginningorendoftheconvolution,sotheresultnowhas20+3-1=22samplesinslowtime.注意此處3脈沖MTI濾波的結(jié)果嚴(yán)格將有誤，僅僅去掉濾波結(jié)果的前后各(s-1)/2個采樣點(diǎn)都不夠。必須在慢時間維的總脈沖個

27、數(shù)20中減去(s-1=2)才是正確結(jié)果！20個脈沖經(jīng)過3脈沖對消后的結(jié)果應(yīng)該為20-3+1=18個脈沖。3脈沖對消實(shí)際上是hn=1-21與慢時間上的各個脈沖進(jìn)行脈沖間濾波，即做內(nèi)積。想想一下，編號為#1,#2,#3的3個脈沖與hn做內(nèi)積后變成了一個脈沖，編號為#2,#3,#4的3個脈沖與hn做內(nèi)積后變成了下一個脈沖.同時應(yīng)注意hn這個滑動窗必'須全部在20個脈沖之內(nèi)才有效，最后#18,#19,#20的脈沖經(jīng)過3脈沖MTI對消后就結(jié)束了。也就是hn的第一個系數(shù)與1#脈沖重合開始算，到hn的最后一個系數(shù)與最后一個脈沖20#脈沖重合時結(jié)束！對于上述ypm(i,:)=conv(h,yp(i,:

28、);的結(jié)果是22個脈沖，要去掉ypm(i,:)前端的(s-1=2)個hn沒有完全沒入的卷積結(jié)果，以及ypm(i,:)后面(s-1=2)個hn沒有完全沒入的卷積結(jié)果，才是最終正確的剩余18個脈沖的處理結(jié)果。Mypm,Nypm=size(ypm);%Dopplerprocessandsquare-lawdetectthewhole%clutter-cancelledarrayanddisplaymeshdisp('putingclutter-cancelledrange-Dopplermap')YPM=fft(conj(ypm').*(hamming(Nypm)*ones(

29、1,Mypm),Lfft);%ypm'是456-by-22的,對慢時間做FFT進(jìn)行相干積累ncy,YPM=git_rotate(YPM.*conj(YPM),Lfft/2);%等效Y=fftshift(YPM.*conj(YPM),1);Ihaveseenmanycaseswhereastudentdoesagoodjobofestimatingfrequeexceptthattheyforgettodoariiftshiftontheirdata,orforgetthattheydiddoanfftshiftontheirdata.YPMdB=db(abs(YPM),'pow

30、er');figureRANGE-DOPPLFtPLOTOFCLUTTEF?-CANCELLEDDATA10Dmesh(doppler,rangep,YPMdB')title('RANGE-DOPPLERPLOTOFCLUTTER-CANCELLEDDATA')ylabel('range(km)')xlabel('velocity(m/s)')levels=(max(YPMdB(:)+-1-5-10-15-20-25);figurecontour(doppler,rangep,YPMdB',levels)title(

31、9;RANGE-DOPPLERCONTOURPLOTOFCLUTTER-CANCELLEDDATA')ylabel('range(km)')xlabel('velocity(m/s)')grid%OK,nowlet'sdoa1Drange-onlyCFARtosetathresholdmap.%DefineCFARwindow23inrange,includingtestcell(cellundertest)and2guardcells(altogether4guardcells),so%18actualaveragingcellscfar=on

32、es(23,1)/18;%我以前不知道CFAR是如何實(shí)現(xiàn)的cfar(10:14)=0;%CellUnderTest單元的index是12,它左、右各有2個guardcells,9個averagingcells%Nowconvolveitinrangedimensionwithrange-DopplermapYPM%togetaveragenoisemap;discardhalf-transientstomaintaindataarraysizeN=zeros(size(YPM);fori=1:Lffttemp=conv(YPM(i,:),cfar);%cfar是23點(diǎn)的h(n),temp是長度

33、為478的行向量456+23-1=478N(i,:)=temp(12:end-11);%用cfar卷積：首先軍出會產(chǎn)生(23-1)/2=11的延遲，即half-transients,將其去掉即從temp(12)開始；其次，輸出的最后(23-1)/2=11個點(diǎn)是部分cfar與信號卷積的結(jié)果，也是half-transients,應(yīng)該去掉。最終剩下的temp結(jié)果是完整的cfar在信號中的結(jié)果，其長度與信號等長，N對噪聲功率的估計(jì)結(jié)果，是256-by-456的與數(shù)據(jù)陣YPM大小相同。end%plottheresultingnoisemapNdB=db(abs(N),'power');f

34、iguremesh(doppler,rangep,NdB')RANGE-DOPPLERPLOTOFCFARNOISEESTIMATEweiocity(ni/s)100從上圖的結(jié)果中可以看出噪聲大約為20dBtitle('RANGE-DOPPLERPLOTOFCFARNOISEESTIMATE')ylabel('range(km),)xlabel('velocity(m/s)')%multiplythenoisebyafactorof16(20dB)andthenthresholdthedata小聲功等的估計(jì)它(7J0)則要求的門跖!可以由估計(jì)到的

35、干擾功率乘以個系效用到.即(7.11)同力T戰(zhàn)動率及相應(yīng)的檢停門限地由油利單兀下圉的臨近單元干擾功等峋平均均用的.明以這種口徜R方法林力為元平均£FAR1匚4CFAR).由卜就動率不是fl!確已知而足通M估計(jì)用到的r所以上式中的尺度內(nèi)內(nèi)口不不產(chǎn)式依位中的尺度因子,其將在7J節(jié)中進(jìn)行推導(dǎo)threshold=16*N;%我以前不知道20dB是從以上估計(jì)出來的detected=NaN*ones(size(YPM);detected(YPM>threshold)=YPM(YPM>threshold);%plottheresultingdetectionmaponadBscalef

36、iguremesh(doppler,rangep,db(detected','power')title('RANGE-DOPPLERDETECTIONMAP')ylabel('range(km)')xlabel('velocity(m/s)')%Tosearchforrangebinswithtargets,firstnoncoherentlyintegrate%acrossthefrequencybinsYPMrange=sum(YPM);%注意YPM是做完相干積累后，取模平方后的結(jié)果！，YPMrange是各個距離單元的

37、噪聲功率%noisefloorestimatebasedonmediantoavoidelevationofestimateby%targetresponsesNrange=median(YPMrange);%YPMrange是各個脈沖經(jīng)過非相干積累的結(jié)果,Nrange是各個距離單元的噪聲功率的中值Trange=8*Nrange;%threshold8x(9DB)abovenoiseestimate%Thisloopidentifieswhichrangebinshavelocalpeaksabovethe%threshold.Italsosetsupavectorforplottingcon

38、veniencetocircle%thepeaksthatarefound.spikesr=;%這一步是尋找距離方向上的峰值marker_r=NaN*ones(1,length(YPMrange);%1-by-456fori=2:Mypm-1%為什么是從YPMrange(2:Mypm-1)?因?yàn)閅PMrange(i)是分別與它左右兩邊的值比較if(YPMrange(i)>YPMrange(i+1)&(YPMrange(i)>YPMrange(i-1)&(YPMrange(i)>Trange)spikesr=spikesr;i,YPMrange(i);%看看如何

39、給spikesr后面按行續(xù)值，這里最終spikesr是3-by-2的marker_r(i)=YPMrange(i);%marker_r的大部分值是NaN,其中有值的是搜索出的峰值的大小endendMspikesr,Nspikesr=size(spikesr);figureplot(rangep,db(YPMrange;Nrange*ones(1,length(YPMrange);.%YPMrange是非相干積累結(jié)果，Nrange是中值Trange*ones(1,length(YPMrange),'power');%Trange是閾值holdonplot(rangep,db(ma

40、rker_r,'power'),'-ro')%marker_r中有值的是搜索出的峰值，marker_r是和YPMrange相同長度的行向量holdoffxlabel('range(km),)ylabel('power(dB),)title('RANGEPEAKS')gridtargets=;%NowfindtheDopplerpeak(s)foreachrangebinhavingatarget(s).Keep%adjoiningDopplervaluesaswelltosupportsubsequentinterpolation

41、.fori=1:Mspikesr%Mspikesr是距離方向搜出的峰值的個數(shù)，是spikesr的行數(shù)，這里Mspikesr=3rb=spikesr(i,1)%currentrangebin注意spikesr是3-by-2的，spikesr的第一列存的是峰值在距離單元的序號%searchinDopplerinthisrangebin.Havetodothismodulothe%DFTsizetoallowforDopplerpeaksnearthespectrumedges.spikesd=口；marker_d=NaN*ones(1,Lfft);%峰值的距離單元上有Lfft=256個多普勒單元，

42、故marker_d是1-by-256的fork=1:Lfftkm1=k-1;%km1是k左邊的那個值ifkm1<1km1=km1+Lfft;%modulotheDFTsizetoallowforDopplerpeaksnearthespectrumedgesendkp1=k+1;%km1是k右邊的那個值ifkp1>Lfftkp1=kp1-Lfft;%modulotheDFTsizetoallowforDopplerpeaksnearthespectrumedgesendif(YPM(k,rb)>YPM(kp1,rb)&(YPM(k,rb)>YPM(km1,rb)

43、&(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值對應(yīng)著多普勒單元上的峰值大小endend%endofloopoverDopplerbinsMspikesd,Nspikesd=size(spikesd);figureplot(doppler,db(YPM(:,rb),Nrange*ones(Lfft,1

44、),'power');%YPM(:,rb)是rb距離單元的峰值對應(yīng)的doppler曲線holdonplot(doppler,db(marker_d,'power'),'-ro')%marker_d是搜索出來的doppler峰值holdoffxlabel('velocity(m/s)')ylabel('power(dB)')title('DOPPLERPEAKSFORRANGEBIN',int2str(rb)grid%ForeachlocalpeakinDoppler,use'peakinte

45、rp'%torefinetheamplitudeandlocationestimateusingaquadratic%interpolationfori=1:Mspikesd%notethatpeakinterpworksonmagnitude,notmagnitude-squaredampdel_k=peakinterp(sqrt(spikesd(i,2:4);targets=targets;amp,rangep(rb),spikesd(i,1)+del_k;endend%endofloopoverrangebinscontainingpeaksMtargets,Ntargets=s

46、ize(targets);fprintf('nnINTERIMPARAMETERSOFDETECTEDTARGETS:n')fprintf('nNumberPowerRange(km)DFTIndex')fori=1:Mtargetsfprintf('n%2.0g%7.3g%6.3g%9.4g',.i,targets(i,:)enddisp('')disp('')%OK,now'targets'containsalist,wehope,ofalltargetswiththe%rangebin,int

47、erpolatedpeakmagnitude,andinterpolatedDopplerbinof%each.Nowweneedtogetinthedesiredunitsandadjusttheamplitudes%foracoupleoffactors.%PutDopplerpeaksintom/sunitsDFTindiceskareconvertedtovelocityusingtherelationshipv=k(/2)(PRF/K).Rememberthatkcountsfrom0inthisformula,whereasMATLABcountsfrom1.targets(:,3

48、)=(targets(:,3)-1)/Lfft)-0.5)*vua;%AdjusttargetamplitudesforMTIfilter,%andthenconverttorelativeRCStargets(:,1)=targets(:,1)./(4*(sin(targets(:,3)*(pi/vua).A2);fprintf('nnINTERIMPARAMETERSOFDETECTEDTARGETSAFTERMTIFILTERCORRECTION:n')fprintf('nNumberPowerRange(km)Vel(m/s)')fori=1:Mtarg

49、etsfprintf('n%2.0g%7.3g%6.3g%9.4g',.i,targets(i,:)enddisp('')disp('')NotethattheR4isthecorrectscalingfortargetdata,butnotforclutter,whichwouldrequireeitherR2orR3,dependingonwhetheritwasbeam-orpulse-limitedclutter.(EarlierInotedthatthecodethatgeneratedthedataassumedbeam-limite

50、dbutyoudon'tknowthatbasedonjusttheinformationprovidedtoyou.)However,IhavefilteredoutmostoftheclutterandIamconcernedwithtargetamplitudes,notclutteramplitude.Also,Iappliedtherangecorrectiontopower(magnitude-squared)data.IfIwantedtodothecorrectiononvoltagedata,IwouldusR2insteadofR4.targets(:,1)=tar

51、gets(:,1).*targets(:,2).A2;targets(:,1)=targets(:,1)/max(targets(:,1);targets(:,1)=db(targets(:,1),'voltage');%Listoutdetectedtargetamplitudes,ranges,Dopplersfprintf('nnESTIMATEDPARAMETERSOFDETECTEDTARGETS:n')fprintf('nNumberRelRCS(dB)Range(km)Vel(m/s)')fori=1:Mtargetsfprintf('n%2.0g%7.3g%6.3g%9.4g',.i,targets(i,:)enddisp('')disp('')用到的函數(shù)functionamp,del_index=peakinterp(z)%PEAKINTERP%peakinterpp