From 29ae9585f97ff18bb9211826f9fd124c74b77e14 Mon Sep 17 00:00:00 2001 From: Ksyer <> Date: Mon, 11 Nov 2024 16:32:53 +0800 Subject: [PATCH] Add Example --- 0 - example/20221025_那斯琦_ICASSP2022.zip | Bin 0 -> 9368 bytes 0 - example/CS-Recovery-Algorithms | 1 + 0 - example/CompressedSensing.jl | 1 + 0 - example/nsq - ICASSP2022/Fig. 2/FISTA.m | 27 ++ ...plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m | 66 +++++ ...test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m | 247 +++++++++++++++++ 0 - example/nsq - ICASSP2022/Fig. 3/FISTA.m | 27 ++ .../plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m | 66 +++++ .../test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m | 249 +++++++++++++++++ 0 - example/nsq - ICASSP2022/Fig. 4/FISTA.m | 27 ++ ...ot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m | 66 +++++ ...st_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m | 251 ++++++++++++++++++ 0 - example/nsq - 恢复算法/FISTA.m | 27 ++ 0 - example/nsq - 恢复算法/cVAMPa_dampling.m | 79 ++++++ .../nsq - 恢复算法/cal_debiased_LASSO.m | 37 +++ .../nsq - 恢复算法/generate_matrix_new.m | 18 ++ 0 - example/nsq - 恢复算法/model.m | 146 ++++++++++ 0 - example/nsq - 恢复算法/sft_thd.m | 20 ++ 0 - example/nsq - 恢复算法/test_model.m | 175 ++++++++++++ 0 - example/wzk - PSK_ROC/test_function.m | 69 +++++ 0 - example/wzk - PSK_ROC/test_signal_t.m | 57 ++++ 0 - example/wzk - PSK_ROC/yAxn_MC.m | 106 ++++++++ 0 - example/wzk - PSK_ROC/yAxn_recovery.m | 57 ++++ 0 - example/wzk - RDA/angle_estimation.m | 60 +++++ 0 - example/wzk - RDA/doppler_process_CS.m | 200 ++++++++++++++ 0 - example/wzk - RDA/doppler_process_MF.m | 45 ++++ 0 - example/wzk - RDA/echo_processing_CS.m | 56 ++++ 0 - example/wzk - RDA/echo_processing_CS_v0.m | 28 ++ 0 - example/wzk - RDA/echo_processing_MF.m | 54 ++++ 0 - example/wzk - RDA/echo_processing_MF_v0.m | 28 ++ 0 - example/wzk - RDA/generate_chirp_mtx.m | 81 ++++++ 0 - example/wzk - RDA/get_range_speed_val.m | 35 +++ 0 - example/wzk - RDA/main_CS.m | 81 ++++++ 0 - example/wzk - RDA/main_MF.m | 79 ++++++ 0 - example/wzk - RDA/plot_result.m | 38 +++ 0 - example/wzk - RDA/range_process_CS.m | 155 +++++++++++ 0 - example/wzk - RDA/range_process_MF.m | 36 +++ 0 - example/wzk - RDA/rda_detection.m | 47 ++++ 0 - example/wzk - RDA/spatial_FFT.m | 31 +++ 0 - example/wzk - main_process/angle_est.m | 43 +++ .../wzk - main_process/doppler_process_CS.m | 187 +++++++++++++ .../wzk - main_process/doppler_process_MF.m | 38 +++ .../wzk - main_process/draw_target_map.m | 21 ++ .../wzk - main_process/generate_chirp_mtx.m | 72 +++++ .../wzk - main_process/get_range_speed_val.m | 29 ++ 0 - example/wzk - main_process/main_cs_part.m | 84 ++++++ 0 - example/wzk - main_process/main_mf.m | 67 +++++ 0 - example/wzk - main_process/main_mf_part.m | 73 +++++ 0 - example/wzk - main_process/plot_result.m | 47 ++++ .../wzk - main_process/range_process_CS.m | 146 ++++++++++ .../wzk - main_process/range_process_MF.m | 30 +++ 0 - example/wzk - main_process/rd_detection.m | 32 +++ .../wzk - 目标检测/generate_chirp_mtx.m | 71 +++++ 0 - example/wzk - 目标检测/node_detect.m | 87 ++++++ 0 - example/wzk - 目标检测/node_detect_norm.m | 87 ++++++ 0 - example/wzk - 目标检测/plot_node_detect.m | 42 +++ 56 files changed, 4029 insertions(+) create mode 100644 0 - example/20221025_那斯琦_ICASSP2022.zip create mode 160000 0 - example/CS-Recovery-Algorithms create mode 160000 0 - example/CompressedSensing.jl create mode 100644 0 - example/nsq - ICASSP2022/Fig. 2/FISTA.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 3/FISTA.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 4/FISTA.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m create mode 100644 0 - example/nsq - ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m create mode 100644 0 - example/nsq - 恢复算法/FISTA.m create mode 100644 0 - example/nsq - 恢复算法/cVAMPa_dampling.m create mode 100644 0 - example/nsq - 恢复算法/cal_debiased_LASSO.m create mode 100755 0 - example/nsq - 恢复算法/generate_matrix_new.m create mode 100755 0 - example/nsq - 恢复算法/model.m create mode 100644 0 - example/nsq - 恢复算法/sft_thd.m create mode 100644 0 - example/nsq - 恢复算法/test_model.m create mode 100644 0 - example/wzk - PSK_ROC/test_function.m create mode 100644 0 - example/wzk - PSK_ROC/test_signal_t.m create mode 100644 0 - example/wzk - PSK_ROC/yAxn_MC.m create mode 100644 0 - example/wzk - PSK_ROC/yAxn_recovery.m create mode 100644 0 - example/wzk - RDA/angle_estimation.m create mode 100644 0 - example/wzk - RDA/doppler_process_CS.m create mode 100644 0 - example/wzk - RDA/doppler_process_MF.m create mode 100644 0 - example/wzk - RDA/echo_processing_CS.m create mode 100644 0 - example/wzk - RDA/echo_processing_CS_v0.m create mode 100644 0 - example/wzk - RDA/echo_processing_MF.m create mode 100644 0 - example/wzk - RDA/echo_processing_MF_v0.m create mode 100644 0 - example/wzk - RDA/generate_chirp_mtx.m create mode 100644 0 - example/wzk - RDA/get_range_speed_val.m create mode 100644 0 - example/wzk - RDA/main_CS.m create mode 100644 0 - example/wzk - RDA/main_MF.m create mode 100644 0 - example/wzk - RDA/plot_result.m create mode 100644 0 - example/wzk - RDA/range_process_CS.m create mode 100644 0 - example/wzk - RDA/range_process_MF.m create mode 100644 0 - example/wzk - RDA/rda_detection.m create mode 100644 0 - example/wzk - RDA/spatial_FFT.m create mode 100644 0 - example/wzk - main_process/angle_est.m create mode 100644 0 - example/wzk - main_process/doppler_process_CS.m create mode 100644 0 - example/wzk - main_process/doppler_process_MF.m create mode 100644 0 - example/wzk - main_process/draw_target_map.m create mode 100644 0 - example/wzk - main_process/generate_chirp_mtx.m create mode 100644 0 - example/wzk - main_process/get_range_speed_val.m create mode 100644 0 - example/wzk - main_process/main_cs_part.m create mode 100644 0 - example/wzk - main_process/main_mf.m create mode 100644 0 - example/wzk - main_process/main_mf_part.m create mode 100644 0 - example/wzk - main_process/plot_result.m create mode 100644 0 - example/wzk - main_process/range_process_CS.m create mode 100644 0 - example/wzk - main_process/range_process_MF.m create mode 100644 0 - example/wzk - main_process/rd_detection.m create mode 100644 0 - example/wzk - 目标检测/generate_chirp_mtx.m create mode 100644 0 - example/wzk - 目标检测/node_detect.m create mode 100644 0 - example/wzk - 目标检测/node_detect_norm.m create mode 100644 0 - example/wzk - 目标检测/plot_node_detect.m diff --git a/0 - example/20221025_那斯琦_ICASSP2022.zip b/0 - example/20221025_那斯琦_ICASSP2022.zip new file mode 100644 index 0000000000000000000000000000000000000000..bbab81aacae3bd11bb3bfd63af1fdc3b552c81d1 GIT binary patch literal 9368 zcmeI2cUTi^xAp@B=@5|KTj;%akkFM55_(5E2t=Al@4Z(ALKgu6=|y@+1Oe#*L5g%~ zQj~-H?05T~^LftZ+TXutawU^Vu377ubx--ty+-3c3Xl*00H6aPa9=13(%X0AH%mZL z0D$pk&d0;Y$IHVfXxh_1+_Z6SstmcOt*!C%S8gsjZBuOz7dPwO3|`DG8Ib&S4=qwQ zC+8V3$R`j4CpkcdBjFZsA?MWqY?c8U}v zhm^}tzFvI0!pj~v@8&L&VvI^PSqI^$?#E5GX=?+nj|PBb(dnc@p9_$+so#~oqB}WX z-b$$qiam^oqY#Wc^-$64>}!pP+X{K|rbN3Lb3NahQYyw@NR4IPR=ory>%oKW2D>)R z^ZY>~d&NV%Q7+N#PHU@qX5eq=ux`|d@WhYk3SRl9Sv4F0Li7@8{8UU>`sjdFd#G+>o6ZIp2>`ghd9&VXeFbH0oqJqPs)?#D!93Vu z`$|`6xnkj60ruj%kWcTGlQp>KonEqubwz!H>ha-x-F}3|kP6x9aWN6cn{v6f`?60h zwL|@4`SZ!HWH)OFM@ao7g;L3qETggnbXS%6`mR@|_rF`n7L!FGpoeBs~}hks+{0C}%WiRFwZA zBZ{17{?B>@=u%SZu zfUcsoqZcazlR1HR4SqZ){DJhV);l+~#_QKXcfG0e8Y%$r;8wMEb#(DC z)vz%))vz+vR?{-oP%wpPsmqx{?x|{+YRf$^{kd@S2ba^7*06my52%&m3a=R`KFO&J4?TRXJAwB7Y>#ayn9xVX)*-l6TzoE3FPXJ7%Kc zXL4GYB|WwBjGJ4u^p7bosRNm31kpT5Xb6Pt?`xp+nUJHZREUfV%E>brdvIXZ%4k00 z8)rbu;6HU4n?m61ezN+W$|5(~4X8qiLQ+G2d%z&+LD&bC-7B|HmkmXt>|b}`Qr}X1 zC&f95d6$?hSap-Dq*0KZKunqsFiQ7IKh)Nhn|(Z~37Q;Q3H8fMfWV3j?FuUc&v}%R zB8)!OcveRXpZe_>JQ9Q8#wiZ_cgA0hJXza-F~6n(f+Hbg6j67voK zU_xj9AzLU$m5a6kmDXJ9iFlvP(bMh2L_N*?tS23!8fj#l__HlbYSg3P$^00<+c=`NzZ70Sl ziEq#ve7`E8zj5@zWBo4VR@(+l!RYgMDPBOb!eBt`9W^k4e3v z>Qb^BRoo**z8WP~ev*fWtt!O|YXh@G>=VG}Xo=7er{lr!+k7u2+PT>B{F6TW4A(6(*o?>(M8Zyobrq-2k3Ld~wv_uzW{)AH z-bF9KYYpMGdhl{>c)2>f<^H7zk*^tBvh0CHk1S&eRu~Lf*VU(M`Q<}oXde_`WTBrd zZQzt?gYgH0TVCEs{RHuB<#3c&>-}t8fA6ODd#OKczIfxfZUU|kMrfts4WOaSIm6XK zsY zB~-_|uE7-)JfpBv2C)u&667d2luQNdh$kffi`+JM_6}psDyoY4*G8qqYfks=0RMYtw z0UP(o8NSuw+}*MW9IR1l(aln>WJ<|$ugd%&)SFFU99}XU++=d9CDBJ4X0CB2HN9FS zDZ@G`C3HG#u49-;*4X6gAucG}A2nxMY82)ubZ|~pnsMIOA(KB#mOWMh*BF7;WchuS z`Wm!g@z$h>P{~vmiOe#p(3&6zFz`d&4+y0_NsAheL@_7zX>X4o^9$YVe(xZ^f7dxw zMNk&cfh9&LK2Al(wd^@vQz;Vjc>&9M`T*xZx*+3qh@q~l-0rZ~M;^OpqdqZIx|&(3 zuzPD102)3**8;6`<1z2#!JghuhXL%Cs(mEl#}qS1FjJeoGy7`xTB)4TDIkQkQ5u`@ zNvYAv!RVb1lPt!PbSf&P_9wvmWo^^-x>{t;bArcNI2?>d?r9H>)&zRiH#I9VO7vV~ z_umM)YO4YbLKer0a9A(svY7fpvFA(h1S%@NvkOGR;^MA-+(P_>5LFwSo0ftmxt>-b zdzs@>?6iu{-(iZIG^;2J6cS@!P7EAf8AS=OdKrK0SesQTOHrD%v5%qwQ*%B0&UY+R zhOAiiX+LwBp z9!ur$d{^`F_Xx~aZ}wS+-B0_#RK6+I#H^18YxDIf-%Dv@PNNjm#wkWjKNf%ycU)6F z+{b0VxDb!wUa;^(!57)7j;n)eE?Uk;tL&mU$6Gf*mibt<)$_-OKwvBa-#Yd};Xe=sn!uT$Huf9@Eop|eOaf*>2I|Ng23fgbqb4TA` zz?kjjl+19|bf%euqhh-C%iKA**%9wJF$L*T;5IcsUsH(A5mFY$8>BBi;sJ5&O%(@Z zB?+FnU>d7FG{DO~PoJkm*L6xis5448%_vFfq{NBg+AMPEN)gMw4K>m*a8`C ziIt5kv8|l7?3z3y({$j_Mv7$3AsizCRXu7A2r=Ij15&$e0CsR`Fon|YzXtX$(HOfG zi13nYJRBsZf`t|5@k?hO@sezcWlj3qPn;fnv0_rOl66OP&j-T8+q$@~RPi8Tr1wkx z^mgP~5argzne;}R26HfNgNCAvXT>|^`S#}DIm};EcYn-sagPXkSA}V&9gsyd;1gd1 z63mPy($5ImLJ22)SrnM8RB-QHAOap>wE7y)J6d?_*`83R(~+3Soc#~e$DcQp*rg%R zbFfvase+DwJK`7J6FWe4Q_Rb=IdAk8q`a-hw_?Db8vnlt2&M|N`V#@$ynZ5}?A*hA z5_K&L=I+yZ&wi4~l+SV|)%mFFw22;1_nKB06D8Oi26v~n2A#SJWQ~rJl3A3OfBa-X zM>xHTcQ2OK*OyVXQ(MX^CIr&oh(TlvsNC)NGrgA%DY?<%K`vcvYN9INLP4Y5){NBs zfnG}wy4kAWK&uaf0*$Y(Ke=|*y7Gt8HR|g+yzGqAP=z0d$;*zD1o4>^X^Jc+=lR<| z7z^JQY6Xndd^xrvDq@`|Gui2N};ZCeqXX2%FxUqPOT_Po?x zp9{iA(RhsXV8m!Qk9lS!hLc!96P+jRRgk*Wow#)JgLxPS%2pIx!n;j+yC$GGW!LrG z{3jU$?#D**gg^-5@f7Uo3Z?=H_5Az@!(MEcfEJzuTZ|Bnm>_Lq&D~_5dl0S@asA^e z!Uh(^WXdtVI&)R`PPXfO0HMjfmStO&Vtua-)Ji5-ItcX?Q*Y}w?<2m<8HsOX1UQ~k zQ{D8=cf~912~JM(Qkn5GnMk`jR(J+pLq11kR|x^uCsKUuX9uZ!bnkX**;m!oM6t!D zwS*bezUT~bs61M$Zg!vYjJPwdiPc!#pk~$>QIAj61-TfqvQotVQZQ+AJpTyKW4vN< z_@iANZrqrp@Reo#68FZkNT0*iRV_+)gL=%VGLq;Qm{!K2mLQw1->my8Nn*$2g_K&6 zAzR}qKut5RzKYfqP_QVTTs^)0;>nXnEhJZ+Et-2XRmUs7=Il5~Th4HUb3m`j7W(M~ zouE9|TD$${-t8JEreUF)bl$Y9e;e19c>) z_V?D(uS}Z6g^#T$rd*qyD4s|1qyN~ih})|O;(pukEj_!uiWd2DQCH8MXmT;NskBGR z%u~d{ngB_a-dO7PvL{6*C>elsB=w=^7$cKak>d4>wsrHfH4B}|h~ z6cQe9iV?on6A2L^bw)noyZw4^=_Za|)pB2Rg=PVR`>WRcjH1;ytt? zQ!BQHj{!1nz!g%W`p%^E0*0gTT!g96X;^5g)yY}L5kWKU7QXA-?f7Ws1viJIzMU`q zwP)hm%~_JqWRWN*O3#s{>^hUjzFjM3fHso9JyyStnk5PLf+xvRxS7oY_OT3a$#cuo z+??7kAFWm88|H*M+NR*Z9b2lZ4v9)H{S`@Ex+m#!!V>u?&};n^P+u)yC3VmHVDhYR ze~toaG%;JYhl^^}V^~DnG_2F&8ar>J?>cnpD4gZc-WLgfFiO_0j<-=|T@h{H7t|X* z>XFpdaB)F9-u^MD8E8mATMh28#d4&YW-@oeorZ(*weO0ZZn*sGs`Sq&QQ-fR5(WN` z6917ucz#2P&5m=M0XLK=b3+kox1dCU+fm{pb0;VBTRMa~e>j9{V}?V8BU?;28XIXqeMyweV#%#IW_@M8sYD8s%H(UYt9jsPjw_K+-3{N$6S)>VWRP`5 zDYg7z(2sNU>7z%=e;S1tNl0qg*RFZgKDd?-=+1Pjne)47KTN&E?^7~P4$~iJLKyKC z(TH9;(lIQflLcei$}_)(f&Z&xXg*-{YEV}1m6XN+UCV@mfh~FlUNB?i0QlKyn&)BP zcKqUW-@CS8%cJ#cTa#z>`(Hmm({sb!&%H^{E-z94bzGtY=x=E8RnERR(~ZuxH?(*= zt#3t(e`@{TfU(9=U!@suRICF!PIBWJ{sD}tvwr|1MaiGQShzN6a086Rw24YmPxrh$ z*Pf@q>f5EH4Ie2R9s24sGrr^Vl9Fy|q`hp6Ix19V(_&UlCZ0d44U(W+DcBdWIJqD; z=mfjy3`*HgNOZnfi}GBQ%gzvfMCYlps!V1Iu4$kkz0~eaUJ-q#BRA>g zE4%Vw#mOO!*8pz+W>Kp6 z7ertc=!uS{r&d2-N%BYQ)c)BX1}t?a&fGB!s`xc(3`#z5r|+k@3|`xuU>iyU=xYzi zflT1{rIQtKH(jGp3!Y|JIOIEBNzf^p+?NpLhs+S2g;uor)iv z<^9~5B=C;=1<`}Ef&4{^dTU(Gy@my+|Ni9N@+`xDdMgTotbUd2R zeqxcLoIMW?{Pmbg|exqxpxM7ukI0_us3rKZE%9Whm>{EfG#DPw`B;p^C7N;;DwR8BZwY;<1$??TjJBj z-3pLd{raeJh+lDdhK-(Ve}j;$3+`F(RRmk05et>5Jb_SCwlRCGK$UWy#H{C9D z_>-6q~_1nKIeFrL?7^<->9!E7Bz+%zCwZHT0~&8je} zERZGma1rc5DXb19CwUX=1RxX4L?ioBMT8ewi; zPvrPXJt23{9)t;%^U$4wZEoUBs>p+J!6{9Jo4Dk#&NOaPY@QV4FAQRougy=!IyL0Q zzq4k09$Lnk5-6cg+eDUHt0*hi&>Mo0lp26b;9b={Fj)X`)GsdG!+L^d;pdl8guCs~ z{lFTkPC~4&<*3?oEnib3Dm`B~kEbhpULw3ReyzL$bD^Nlw6aTFD0(%k$~aazLdse# zxA?H^sh9$i11VetIUF@shDMY7emw}&Bj^hR0fcpkCyZqkIv delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = mean(abs(z_pre - z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; +end + +end \ No newline at end of file diff --git a/0 - example/nsq - ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m b/0 - example/nsq - ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m new file mode 100644 index 0000000..dc7b77f --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m @@ -0,0 +1,66 @@ +clear; +close all; +clc; + +load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.mat; + +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + +%% plot +figure(1); +plot(SNR, P_fa_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(SNR, P_fa_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(SNR, P_fa_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(SNR, P_fa_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('SNR'); +ylabel('P_{fa}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + +figure(2); +plot(SNR, P_d_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(SNR, P_d_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(SNR, P_d_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(SNR, P_d_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('SNR'); +ylabel('P_{d}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + + + + + + + + + + diff --git a/0 - example/nsq - ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m b/0 - example/nsq - ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m new file mode 100644 index 0000000..e018e4d --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m @@ -0,0 +1,247 @@ +clc; +clear; +close all; + +%% parameter setting +m = 128; +n = 256; + +SNR = 0: 1: 15; +len_SNR = length(SNR); + +rep_time = 1e4; + +P_fa = 1e-2; + +p0 = 0.1; + +lambda = 0.1; + +sigma_0 = 0.1; + +sigma_n = sigma_0; + +%% experiment +gamma = m/n; + +P_fa_CROD_cnt = zeros(len_SNR, rep_time); +P_fa_CAMP_cnt = zeros(len_SNR, rep_time); +P_fa_SDL_cnt = zeros(len_SNR, rep_time); +P_fa_ROD_cnt = zeros(len_SNR, rep_time); +P_fa_LASSO_cnt = zeros(len_SNR, rep_time); + +P_d_CROD_cnt = zeros(len_SNR, rep_time); +P_d_CAMP_cnt = zeros(len_SNR, rep_time); +P_d_SDL_cnt = zeros(len_SNR, rep_time); +P_d_ROD_cnt = zeros(len_SNR, rep_time); +P_d_LASSO_cnt = zeros(len_SNR, rep_time); + +x_idx = rand(n, 1); +if p0 == 0 + thd = -1; + x_l0 = sum(x_idx > thd); +else + thd = sort(x_idx); + thd = thd(round(n*p0)); + x_l1 = sum(x_idx <= thd); + x_l0 = sum(x_idx > thd); +end + +x = zeros(n, 1); +x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1)); +x(x_idx <= thd) = x_temp(x_idx <= thd); +x = x * sqrt(n/m); + +h_thd = -log(P_fa); + +parfor rep = 1: rep_time + + A_idx = randperm(n); + A_idx = A_idx(1: m); + A_idx = sort(A_idx); + A = dftmtx(n); + A = A(A_idx, :); + A = A / sqrt(n); + + w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1); +% w = w / sqrt(10^(SNR(cnt_SNR)/10)); + + for cnt_SNR = 1: len_SNR + + x1 = x * sqrt(10^(SNR(cnt_SNR)/10)); + + y = A * x1 + w; + + % x_LASSO = LASSO_cvx(y, A, lambda); + x_LASSO = FISTA(y, A, lambda, 1e-5); + + % CROD + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat; + RSS = sum(abs(y - A * x_LASSO).^2)/m; + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + stat_CROD = abs(x_d_CROD / sigma_CROD).^2; + + P_fa_CROD_cnt(cnt_SNR, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0; + P_d_CROD_cnt(cnt_SNR, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1; + + % CAMP + Q_hat1 = gamma - rho_active; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat1 = (gamma-Rho); + x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1; + sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP)); + stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2; + + P_fa_CAMP_cnt(cnt_SNR, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0; + P_d_CAMP_cnt(cnt_SNR, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1; + + % SDL + Q_hat2 = (gamma - rho_active); + x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2; + sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO)); + stat_SDL = abs(x_d_SDL / sigma_SDL).^2; + + P_fa_SDL_cnt(cnt_SNR, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0; + P_d_SDL_cnt(cnt_SNR, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1; + + % ROD + Q_hat3 = (gamma - rho_active)/(1 - rho_active); + x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3; + chi = rho_active*(1 - rho_active)/(gamma - rho_active); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_ROD = sqrt(2*chi_hat2) / Q_hat3; + stat_ROD = abs(x_d_ROD / sigma_ROD).^2; + + P_fa_ROD_cnt(cnt_SNR, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0; + P_d_ROD_cnt(cnt_SNR, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1; + + % LASSO +% stat_LASSO = abs(x_LASSO).^2; +% if p0 ~= 0 +% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd); +% end +% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd); + + + end + + fprintf('%d\n', rep); + +end + +P_fa_CROD = mean(P_fa_CROD_cnt, 2); +P_fa_CAMP = mean(P_fa_CAMP_cnt, 2); +P_fa_SDL = mean(P_fa_SDL_cnt, 2); +P_fa_ROD = mean(P_fa_ROD_cnt, 2); + +P_d_CROD = mean(P_d_CROD_cnt, 2); +P_d_CAMP = mean(P_d_CAMP_cnt, 2); +P_d_SDL = mean(P_d_SDL_cnt, 2); +P_d_ROD = mean(P_d_ROD_cnt, 2); + + +%% plot +figure(1); +plot(SNR, P_fa_CROD, 'linewidth', 2); +hold on; +grid on; +plot(SNR, P_fa_CAMP, 'linewidth', 2); +plot(SNR, P_fa_SDL, 'linewidth', 2); +plot(SNR, P_fa_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('SNR'); +ylabel('P_{fa}'); + + +figure(2); +plot(SNR, P_d_CROD, 'linewidth', 2); +hold on; +grid on; +plot(SNR, P_d_CAMP, 'linewidth', 2); +plot(SNR, P_d_SDL, 'linewidth', 2); +plot(SNR, P_d_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('SNR'); +ylabel('P_{d}'); + + +save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.mat ... + SNR... + P_fa_CROD... + P_fa_CAMP... + P_fa_SDL... + P_fa_ROD... + P_d_CROD... + P_d_CAMP... + P_d_SDL... + P_d_ROD; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - ICASSP2022/Fig. 3/FISTA.m b/0 - example/nsq - ICASSP2022/Fig. 3/FISTA.m new file mode 100644 index 0000000..24f7223 --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 3/FISTA.m @@ -0,0 +1,27 @@ +function [z] = FISTA(y, A, lambda, delta) + +x_pre = A'*y; +t = 1; +z = x_pre; +z_pre = z; +t_pre = t; +N = size(A, 2); +diff = 1; +E = eig(A'*A); +L = E(end); +temp1 = A'*y/L; +temp2 = eye(N) - A'*A/L; +k = 0; +while((diff > delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = mean(abs(z_pre - z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; +end + +end \ No newline at end of file diff --git a/0 - example/nsq - ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m b/0 - example/nsq - ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m new file mode 100644 index 0000000..40ef17b --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m @@ -0,0 +1,66 @@ +clear; +close all; +clc; + +load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.mat; + +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + +%% plot +figure(1); +plot(p0_total, P_fa_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(p0_total, P_fa_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(p0_total, P_fa_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(p0_total, P_fa_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('signal density'); +ylabel('P_{fa}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + +figure(2); +plot(p0_total, P_d_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(p0_total, P_d_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(p0_total, P_d_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(p0_total, P_d_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('signal density'); +ylabel('P_{d}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + + + + + + + + + + diff --git a/0 - example/nsq - ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m b/0 - example/nsq - ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m new file mode 100644 index 0000000..b436281 --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m @@ -0,0 +1,249 @@ +clc; +clear; +close all; + +%% parameter setting +m = 128; +n = 256; + +SNR = 13; + +rep_time = 1e4; + +P_fa = 1e-2; + +p0_total = 0.02: 0.02: 0.2; +len_p0 = length(p0_total); + +lambda = 0.1; + +sigma_0 = 0.1; + +sigma_n = sigma_0; + +%% experiment +gamma = m/n; + +P_fa_CROD_cnt = zeros(len_p0, rep_time); +P_fa_CAMP_cnt = zeros(len_p0, rep_time); +P_fa_SDL_cnt = zeros(len_p0, rep_time); +P_fa_ROD_cnt = zeros(len_p0, rep_time); +P_fa_LASSO_cnt = zeros(len_p0, rep_time); + +P_d_CROD_cnt = zeros(len_p0, rep_time); +P_d_CAMP_cnt = zeros(len_p0, rep_time); +P_d_SDL_cnt = zeros(len_p0, rep_time); +P_d_ROD_cnt = zeros(len_p0, rep_time); +P_d_LASSO_cnt = zeros(len_p0, rep_time); + +h_thd = -log(P_fa); + +parfor rep = 1: rep_time + + A_idx = randperm(n); + A_idx = A_idx(1: m); + A_idx = sort(A_idx); + A = dftmtx(n); + A = A(A_idx, :); + A = A / sqrt(n); + + w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1); +% w = w / sqrt(10^(SNR(cnt_SNR)/10)); + + for cnt_p0 = 1: len_p0 + + p0 = p0_total(cnt_p0); + + x_idx = rand(n, 1); + if p0 == 0 + thd = -1; + x_l0 = sum(x_idx > thd); + else + thd = sort(x_idx); + thd = thd(round(n*p0)); + x_l1 = sum(x_idx <= thd); + x_l0 = sum(x_idx > thd); + end + + x = zeros(n, 1); + x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1)); + x(x_idx <= thd) = x_temp(x_idx <= thd); + x = x * sqrt(n/m); + + x1 = x * sqrt(10^(SNR/10)); + + y = A * x1 + w; + + % x_LASSO = LASSO_cvx(y, A, lambda); + x_LASSO = FISTA(y, A, lambda, 1e-5); + + % CROD + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat; + RSS = sum(abs(y - A * x_LASSO).^2)/m; + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + stat_CROD = abs(x_d_CROD / sigma_CROD).^2; + + P_fa_CROD_cnt(cnt_p0, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0; + P_d_CROD_cnt(cnt_p0, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1; + + % CAMP + Q_hat1 = gamma - rho_active; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat1 = (gamma-Rho); + x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1; + sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP)); + stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2; + + P_fa_CAMP_cnt(cnt_p0, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0; + P_d_CAMP_cnt(cnt_p0, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1; + + % SDL + Q_hat2 = (gamma - rho_active); + x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2; + sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO)); + stat_SDL = abs(x_d_SDL / sigma_SDL).^2; + + P_fa_SDL_cnt(cnt_p0, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0; + P_d_SDL_cnt(cnt_p0, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1; + + % ROD + Q_hat3 = (gamma - rho_active)/(1 - rho_active); + x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3; + chi = rho_active*(1 - rho_active)/(gamma - rho_active); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_ROD = sqrt(2*chi_hat2) / Q_hat3; + stat_ROD = abs(x_d_ROD / sigma_ROD).^2; + + P_fa_ROD_cnt(cnt_p0, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0; + P_d_ROD_cnt(cnt_p0, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1; + + % LASSO +% stat_LASSO = abs(x_LASSO).^2; +% if p0 ~= 0 +% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd); +% end +% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd); + + + end + + fprintf('%d\n', rep); + +end + +P_fa_CROD = mean(P_fa_CROD_cnt, 2); +P_fa_CAMP = mean(P_fa_CAMP_cnt, 2); +P_fa_SDL = mean(P_fa_SDL_cnt, 2); +P_fa_ROD = mean(P_fa_ROD_cnt, 2); + +P_d_CROD = mean(P_d_CROD_cnt, 2); +P_d_CAMP = mean(P_d_CAMP_cnt, 2); +P_d_SDL = mean(P_d_SDL_cnt, 2); +P_d_ROD = mean(P_d_ROD_cnt, 2); + + +%% plot +figure(1); +plot(p0_total, P_fa_CROD, 'linewidth', 2); +hold on; +grid on; +plot(p0_total, P_fa_CAMP, 'linewidth', 2); +plot(p0_total, P_fa_SDL, 'linewidth', 2); +plot(p0_total, P_fa_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('signal density'); +ylabel('P_{fa}'); + + +figure(2); +plot(p0_total, P_d_CROD, 'linewidth', 2); +hold on; +grid on; +plot(p0_total, P_d_CAMP, 'linewidth', 2); +plot(p0_total, P_d_SDL, 'linewidth', 2); +plot(p0_total, P_d_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('signal density'); +ylabel('P_{d}'); + + +save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.mat ... + p0_total... + P_fa_CROD... + P_fa_CAMP... + P_fa_SDL... + P_fa_ROD... + P_d_CROD... + P_d_CAMP... + P_d_SDL... + P_d_ROD; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - ICASSP2022/Fig. 4/FISTA.m b/0 - example/nsq - ICASSP2022/Fig. 4/FISTA.m new file mode 100644 index 0000000..24f7223 --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 4/FISTA.m @@ -0,0 +1,27 @@ +function [z] = FISTA(y, A, lambda, delta) + +x_pre = A'*y; +t = 1; +z = x_pre; +z_pre = z; +t_pre = t; +N = size(A, 2); +diff = 1; +E = eig(A'*A); +L = E(end); +temp1 = A'*y/L; +temp2 = eye(N) - A'*A/L; +k = 0; +while((diff > delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = mean(abs(z_pre - z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; +end + +end \ No newline at end of file diff --git a/0 - example/nsq - ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m b/0 - example/nsq - ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m new file mode 100644 index 0000000..baf4bff --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m @@ -0,0 +1,66 @@ +clear; +close all; +clc; + +load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.mat; + +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + +%% plot +figure(1); +plot(gamma_total, P_fa_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(gamma_total, P_fa_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(gamma_total, P_fa_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(gamma_total, P_fa_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('compression rate'); +ylabel('P_{fa}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + +figure(2); +plot(gamma_total, P_d_CROD, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +hold on; +grid on; +plot(gamma_total, P_d_CAMP, '-d', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(gamma_total, P_d_SDL, '-s', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +plot(gamma_total, P_d_ROD, '-+', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('compression rate'); +ylabel('P_{d}'); +set(gca, 'FontSize', Fontsize); +%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman'); +set(gcf, 'position', [200, 300, plot_width, plot_height]); + + + + + + + + + + diff --git a/0 - example/nsq - ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m b/0 - example/nsq - ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m new file mode 100644 index 0000000..212a723 --- /dev/null +++ b/0 - example/nsq - ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m @@ -0,0 +1,251 @@ +clc; +clear; +close all; + +%% parameter setting +n = 256; + +SNR = 13; + +rep_time = 1e4; + +P_fa = 1e-2; + +p0 = 0.1; + +gamma_total = (4: 12)/16; +len_gamma = length(gamma_total); + +lambda = 0.1; + +sigma_0 = 0.1; + +sigma_n = sigma_0; + +%% experiment + +P_fa_CROD_cnt = zeros(len_gamma, rep_time); +P_fa_CAMP_cnt = zeros(len_gamma, rep_time); +P_fa_SDL_cnt = zeros(len_gamma, rep_time); +P_fa_ROD_cnt = zeros(len_gamma, rep_time); +P_fa_LASSO_cnt = zeros(len_gamma, rep_time); + +P_d_CROD_cnt = zeros(len_gamma, rep_time); +P_d_CAMP_cnt = zeros(len_gamma, rep_time); +P_d_SDL_cnt = zeros(len_gamma, rep_time); +P_d_ROD_cnt = zeros(len_gamma, rep_time); +P_d_LASSO_cnt = zeros(len_gamma, rep_time); + +h_thd = -log(P_fa); + +parfor rep = 1: rep_time + + + + for cnt_gamma = 1: len_gamma + + gamma = gamma_total(cnt_gamma); + m = round(gamma*n); + + A_idx = randperm(n); + A_idx = A_idx(1: m); + A_idx = sort(A_idx); + A = dftmtx(n); + A = A(A_idx, :); + A = A / sqrt(n); + + w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1); + + x_idx = rand(n, 1); + if p0 == 0 + thd = -1; + x_l0 = sum(x_idx > thd); + else + thd = sort(x_idx); + thd = thd(round(n*p0)); + x_l1 = sum(x_idx <= thd); + x_l0 = sum(x_idx > thd); + end + + x = zeros(n, 1); + x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1)); + x(x_idx <= thd) = x_temp(x_idx <= thd); + x = x * sqrt(n/m); + + x1 = x * sqrt(10^(SNR/10)); + + y = A * x1 + w; + + % x_LASSO = LASSO_cvx(y, A, lambda); + x_LASSO = FISTA(y, A, lambda, 1e-5); + + % CROD + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat; + RSS = sum(abs(y - A * x_LASSO).^2)/m; + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + stat_CROD = abs(x_d_CROD / sigma_CROD).^2; + + P_fa_CROD_cnt(cnt_gamma, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0; + P_d_CROD_cnt(cnt_gamma, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1; + + % CAMP + Q_hat1 = gamma - rho_active; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat1 = (gamma-Rho); + x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1; + sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP)); + stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2; + + P_fa_CAMP_cnt(cnt_gamma, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0; + P_d_CAMP_cnt(cnt_gamma, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1; + + % SDL + Q_hat2 = (gamma - rho_active); + x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2; + sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO)); + stat_SDL = abs(x_d_SDL / sigma_SDL).^2; + + P_fa_SDL_cnt(cnt_gamma, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0; + P_d_SDL_cnt(cnt_gamma, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1; + + % ROD + Q_hat3 = (gamma - rho_active)/(1 - rho_active); + x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3; + chi = rho_active*(1 - rho_active)/(gamma - rho_active); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_ROD = sqrt(2*chi_hat2) / Q_hat3; + stat_ROD = abs(x_d_ROD / sigma_ROD).^2; + + P_fa_ROD_cnt(cnt_gamma, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0; + P_d_ROD_cnt(cnt_gamma, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1; + + % LASSO +% stat_LASSO = abs(x_LASSO).^2; +% if p0 ~= 0 +% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd); +% end +% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd); + + + end + + fprintf('%d\n', rep); + +end + +P_fa_CROD = mean(P_fa_CROD_cnt, 2); +P_fa_CAMP = mean(P_fa_CAMP_cnt, 2); +P_fa_SDL = mean(P_fa_SDL_cnt, 2); +P_fa_ROD = mean(P_fa_ROD_cnt, 2); + +P_d_CROD = mean(P_d_CROD_cnt, 2); +P_d_CAMP = mean(P_d_CAMP_cnt, 2); +P_d_SDL = mean(P_d_SDL_cnt, 2); +P_d_ROD = mean(P_d_ROD_cnt, 2); + + +%% plot +figure(1); +plot(gamma_total, P_fa_CROD, 'linewidth', 2); +hold on; +grid on; +plot(gamma_total, P_fa_CAMP, 'linewidth', 2); +plot(gamma_total, P_fa_SDL, 'linewidth', 2); +plot(gamma_total, P_fa_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('compression rate'); +ylabel('P_{fa}'); + + +figure(2); +plot(gamma_total, P_d_CROD, 'linewidth', 2); +hold on; +grid on; +plot(gamma_total, P_d_CAMP, 'linewidth', 2); +plot(gamma_total, P_d_SDL, 'linewidth', 2); +plot(gamma_total, P_d_ROD, 'linewidth', 2); +legend('CROD', 'CAMP', 'SDL-test', 'ROD'); +xlabel('compression rate'); +ylabel('P_{d}'); + + +save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.mat ... + gamma_total... + P_fa_CROD... + P_fa_CAMP... + P_fa_SDL... + P_fa_ROD... + P_d_CROD... + P_d_CAMP... + P_d_SDL... + P_d_ROD; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恢复算法/FISTA.m b/0 - example/nsq - 恢复算法/FISTA.m new file mode 100644 index 0000000..736bcd1 --- /dev/null +++ b/0 - example/nsq - 恢复算法/FISTA.m @@ -0,0 +1,27 @@ +function [z] = FISTA(y, A, lambda, delta) + +x_pre = A'*y; +t = 1; +z = x_pre; +z_pre = z; +t_pre = t; +N = size(A, 2); +diff = 1; +E = eig(A'*A); +L = E(end); +temp1 = A'*y/L; +temp2 = eye(N) - A'*A/L; +k = 0; +while((diff > delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = sum(abs(z_pre - z))/sum(abs(z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; +end + +end \ No newline at end of file diff --git a/0 - example/nsq - 恢复算法/cVAMPa_dampling.m b/0 - example/nsq - 恢复算法/cVAMPa_dampling.m new file mode 100644 index 0000000..3194047 --- /dev/null +++ b/0 - example/nsq - 恢复算法/cVAMPa_dampling.m @@ -0,0 +1,79 @@ +function[x, x_d, hat_Q1, sigma_d, ifcvg] = cVAMPa_dampling(y, A, lambda, alpha, delta, iter_max, sigma) + +[M, N] = size(A); + +gamma = M/N; +p = A'*y; +h1 = p/gamma; +hat_Q1 = gamma; + +%Eigenvalue Decomposition +[V, D] = eig(A'*A); +d = diag(D); + +t = 0; +diff = 1; + +while((diff > delta) && (t < iter_max)) + + h1_pre = h1; + + % Factorized + x1 = sft_thd(h1, lambda/hat_Q1); + chi1 = sum((abs(x1) > 1e-4).* (2 - lambda./((hat_Q1*abs(x1) + lambda)))) / 2 / N / hat_Q1; + % Message F to G + hat_Q2 = 1/chi1 - hat_Q1; + h2 = (x1/chi1 - h1*hat_Q1)/hat_Q2; + % Gaussian + tmp = V'*(p + h2*hat_Q2); + tmp = tmp./(d+hat_Q2); + x2 = V*tmp; + chi2 = sum(1./(d+hat_Q2))/N; + % Message G to F + hat_Q1 = 1/chi2 - hat_Q2; + h1 = alpha*(x2/chi2 - h2*hat_Q2)/hat_Q1+(1-alpha)*h1; + + diff = sum(abs(h1_pre - h1))/sum(abs(h1)); + t = t+1; + +end + +ifcvg = diff <= delta; + +x = x1; +x_d = h1; + +chi = chi1; +t = -hat_Q2; +t_prime = -1/mean((1./(d+hat_Q2)).^2); +G_prime = t + 1/chi; +G_wprime = t_prime + 1/chi/chi; +RSS = sum(abs(y - A*x).^2)/M; +hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +... + (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2; +sigma_d = sqrt(2*hat_chi)/hat_Q1; + +end + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恢复算法/cal_debiased_LASSO.m b/0 - example/nsq - 恢复算法/cal_debiased_LASSO.m new file mode 100644 index 0000000..88ed7c6 --- /dev/null +++ b/0 - example/nsq - 恢复算法/cal_debiased_LASSO.m @@ -0,0 +1,37 @@ +function [x_d, hat_Q1, sigma_d] = cal_debiased_LASSO(x, A, y, lambda, sigma) + +[M, N] = size(A); +gamma = M/N; +hat_Q1 = gamma; +[~, D] = eig(A'*A); +d = diag(D); + +diff = 1; +T = 1000; +t = 0; + +while (t < T) && (diff > 1e-6) + + Q1_pre = hat_Q1; + rho = mean((2 - lambda./(hat_Q1*abs(x) + lambda)).*(abs(x) > 1e-4))/2; + hat_Q1 = rho/mean(1./(d + (1-rho)*hat_Q1/rho)); + diff = abs(Q1_pre - hat_Q1); + t = t+1; + +end + +x_d = x + 1/hat_Q1*A'*(y - A*x); + +chi = rho/hat_Q1; +hat_Q2 = 1/chi - hat_Q1; + +t = -hat_Q2; +t_prime = -1/mean((1./(d+hat_Q2)).^2); +G_prime = t + 1/chi; +G_wprime = t_prime + 1/chi/chi; +RSS = sum(abs(y - A*x).^2)/M; +hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +... + (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2; +sigma_d = sqrt(2*hat_chi)/hat_Q1; + +end \ No newline at end of file diff --git a/0 - example/nsq - 恢复算法/generate_matrix_new.m b/0 - example/nsq - 恢复算法/generate_matrix_new.m new file mode 100755 index 0000000..cb47ff7 --- /dev/null +++ b/0 - example/nsq - 恢复算法/generate_matrix_new.m @@ -0,0 +1,18 @@ +function [ A ] = generate_matrix_new(y1, y2) + A = []; + l = length(y1); + temp1 = y1; + temp2 = y2; + for i = 1:l + t1 = circshift(temp1, i-1); + t2 = circshift(temp2, i-1); + if i - 1 > 0 + t1(1:i - 1,1) = 0; + end + if i - 1 > 0 + t2(1:i - 1,1) = 0; + end + A = [A,t1,t2]; + end +end + diff --git a/0 - example/nsq - 恢复算法/model.m b/0 - example/nsq - 恢复算法/model.m new file mode 100755 index 0000000..f56685f --- /dev/null +++ b/0 - example/nsq - 恢复算法/model.m @@ -0,0 +1,146 @@ +clc; +clear; +close all; + +%% 参数设置 +sigma_n = 0.1; +gamma = 0.5; + +% 信号参数 +B = 5e5; %信号带宽 +Tp = 100e-6; %脉宽100us +fs = 2 * B; %采样频率 +Ts = 1 / fs; %采样周期 +K = B / Tp; %线性调频率 +fc = 1e8; %载波频率 + +Tr = 1e-3; + +t = 0: 1/fs: Tr - 1/fs; +t2 = 0: 1/fs/2: Tr - 1/fs/2; + +c = 3e8; % 光速 + +distance_max = (Tr-Tp) * c / 2; + +target_scattering = [0.8, 1, 0.9]; %扩展目标各点散射强度 + + +%% 生成矩阵 A +% 生成发射信号 signal_t 及 +N = Tr * fs; +N_high = Tp * fs; + +signal_t = zeros(1, N); +signal_td = zeros(1, N); +for i = 1:N_high + tp = (i - 1) * (1 / fs); + signal_t(1, i) = exp(1j*2*pi*(fc*tp+0.5*K*tp.^2)); + tp2 = (i - 0.5) * (1 / fs); + signal_td(1, i + 1) = exp(1j*2*pi*(fc*tp2+0.5*K*tp2.^2)); +end + +A = generate_matrix_new(signal_t.', signal_td.'); + +J = A'*A; +lambda_J=eig(J); + +figure(1) +spy(A) + + +%% 生成回波 y +% 设置目标 - 扩展目标，由三个点组成 +distance1 = 52000; +tau1 = distance1 * 2 / c; +n_tau1 = round(tau1 * fs); + +alpha1 = 0.6; % 扩展目标整体散射强度 +signal_r1_t1 = zeros(1, N); +signal_r2_t1 = zeros(1, N); +signal_r3_t1 = zeros(1, N); + +for i = 1:N + temp = i - n_tau1; + if temp >= 1 && temp <= N_high + signal_r1_t1(1, i) = alpha1 * target_scattering(1) * signal_t(1, temp); + if i + 1 <= N + signal_r2_t1(1, i + 1) = alpha1 * target_scattering(2) * signal_t(1, temp); + end + if i + 1 <= N + signal_r3_t1(1, i + 2) = alpha1 * target_scattering(3) * signal_t(1, temp); + end + end +end +signal_r1 = signal_r1_t1 + signal_r2_t1 + signal_r3_t1; + +% 回波 +signal_r = signal_r1; + +% 加入噪声 +noise = random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)); +signal_r_n = signal_r + noise; + +y = signal_r_n.'; +figure(2) +subplot(211); +plot(t, real(signal_t)); +title('发射信号') +xlabel('时间') +ylabel('幅度') + +subplot(212); +plot(t, real(y)) +title('接收信号(y)') +xlabel('时间') +ylabel('幅度') + +%% 理论 x +distance_node = round((distance1 * 2 / c) * fs); +x_t = zeros(1, 2 * N); +for i = 1: length(target_scattering) + x_t((distance_node + i) * 2 - 1) = alpha1 * target_scattering(i); +end +x = x_t.'; + +figure(3) +plot(t2, x); +title('目标散射点(x)') +xlabel('时间') +ylabel('幅度') + + + +%% 验证 +y_t = A * x; +y_r = signal_r.'; + +figure(4) +subplot(311); +plot(real(y_r)) +title('实际回波') + +subplot(312); +plot(real(y_t)); +title('计算结果') + +subplot(313); +plot(abs(y_r - y_t)); +title('差异') + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恢复算法/sft_thd.m b/0 - example/nsq - 恢复算法/sft_thd.m new file mode 100644 index 0000000..95e46ba --- /dev/null +++ b/0 - example/nsq - 恢复算法/sft_thd.m @@ -0,0 +1,20 @@ +function y = sft_thd(x, thd) + +if isequal(size(x), size(thd)) + + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd); + +else + + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd); + +end + + +end \ No newline at end of file diff --git a/0 - example/nsq - 恢复算法/test_model.m b/0 - example/nsq - 恢复算法/test_model.m new file mode 100644 index 0000000..bc23a0e --- /dev/null +++ b/0 - example/nsq - 恢复算法/test_model.m @@ -0,0 +1,175 @@ +clc; +clear; +%close all; + +%% 参数设置 +sigma_n = 0.1; +gamma = 0.5; + +% 信号参数 +B = 5e5; %信号带宽 +Tp = 100e-6; %脉宽100us +fs = 2 * B; %采样频率 +Ts = 1 / fs; %采样周期 +K = B / Tp; %线性调频率 +fc = 1e8; %载波频率 + +Tr = 1e-3; + +t = 0: 1/fs: Tr - 1/fs; +t2 = 0: 1/fs/2: Tr - 1/fs/2; + +c = 3e8; % 光速 + +distance_max = (Tr-Tp) * c / 2; + +target_scattering = [0.8, 1, 0.9]; %扩展目标各点散射强度 + + +%% 生成矩阵 A +% 生成发射信号 signal_t 及 +N = Tr * fs; +N_high = Tp * fs; + +signal_t = zeros(1, N); +signal_td = zeros(1, N); +for i = 1:N_high + tp = (i - 1) * (1 / fs); + signal_t(1, i) = exp(1j*2*pi*(fc*tp+0.5*K*tp.^2)); + tp2 = (i - 0.5) * (1 / fs); + signal_td(1, i + 1) = exp(1j*2*pi*(fc*tp2+0.5*K*tp2.^2)); +end + +A = generate_matrix_new(signal_t.', signal_td.'); + +temp = 0; +for i = 1: size(A, 1) + + for j = 1: size(A, 2) + + temp = temp + abs(A(mod(i, size(A, 1))+1, mod(j+1, size(A, 2))+1) - A(i, j)); + + end + +end + +%% 生成回波 y +% 设置目标 - 扩展目标，由三个点组成 +distance1 = 52000; +tau1 = distance1 * 2 / c; +n_tau1 = round(tau1 * fs); + +alpha1 = 0.6; % 扩展目标整体散射强度 +signal_r1_t1 = zeros(1, N); +signal_r2_t1 = zeros(1, N); +signal_r3_t1 = zeros(1, N); + +for i = 1:N + temp = i - n_tau1; + if temp >= 1 && temp <= N_high + signal_r1_t1(1, i) = alpha1 * target_scattering(1) * signal_t(1, temp); + if i + 1 <= N + signal_r2_t1(1, i + 1) = alpha1 * target_scattering(2) * signal_t(1, temp); + end + if i + 1 <= N + signal_r3_t1(1, i + 2) = alpha1 * target_scattering(3) * signal_t(1, temp); + end + end +end +signal_r1 = signal_r1_t1 + signal_r2_t1 + signal_r3_t1; + +% 回波 +signal_r = signal_r1; + +% 加入噪声 +noise = random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)); +signal_r_n = signal_r + noise; + +y = signal_r_n.'; + +%% 理论 x +distance_node = round((distance1 * 2 / c) * fs); +x_t = zeros(1, 2 * N); +for i = 1: length(target_scattering) + x_t((distance_node + i) * 2 - 1) = alpha1 * target_scattering(i); +end +x = x_t.'; + +%% Experiment + +%% Parameters setting +lambda = 0.002; +alpha = 1/4; +delta = 1e-8*alpha; +iter_max = round(2000/alpha); +n = size(A, 2); + +% J = A'*A; +J1 = A*A'; +lambda_J=eig(J1); +% histogram(lambda_J, 100); + +%% Normalized +A = A / sqrt(lambda_J(end)); +y = y / sqrt(lambda_J(end)); +sigma_n = sigma_n / sqrt(lambda_J(end)); + +%% cVAMP +tic; +[x_VAMP, x_d, hat_Q1, sigma_d, ifcvg] = cVAMPa_dampling(y, A, lambda, alpha, delta, iter_max, sigma_n); +toc; + +%% cvx +tic; +cvx_begin quiet + variable x_cvx(n, 1) complex + z = lambda*sum(abs(x_cvx)) + 0.5*sum(pow_abs((y - A * x_cvx), 2)); + minimize(z) +cvx_end + +[x_d_cal, hat_Q1_cal, sigma_d_cal] = cal_debiased_LASSO(x_cvx, A, y, lambda, sigma_n); +toc; + +sigma_ex = std(x_d_cal - x, 1); +tmp = (x_d_cal - x)/sigma_ex; + +[h_r, p_r, k_r, c_r] = kstest(real(tmp)*sqrt(2)); +[h_i, p_i, k_i, c_i] = kstest(imag(tmp)*sqrt(2)); + +%% results +% whether cVAMP algorithm converges +ifcvg + +% whether the output of cVAMP converges to the LASSO solution +sum(abs(x_cvx - x_VAMP)) + +% check the results from cVAMP and "calculation" +abs(hat_Q1 - hat_Q1_cal) +abs(sigma_d - sigma_d_cal) +sum(abs(x_d - x_d_cal)) + +% accuracy of estimating the variance +abs(sigma_d_cal - sigma_ex)/abs(sigma_ex) + +% p-value of KS-test +% the larger, the higher probability it is drawn from Gaussian distribution +p_r +p_i + + + + + + + + + + + + + + + + + + diff --git a/0 - example/wzk - PSK_ROC/test_function.m b/0 - example/wzk - PSK_ROC/test_function.m new file mode 100644 index 0000000..dd625be --- /dev/null +++ b/0 - example/wzk - PSK_ROC/test_function.m @@ -0,0 +1,69 @@ + +% load st.mat +% [ A ] = generate_matrix_long(transpose(s_T)); + + +N = 512; +A = dftmtx(N) / sqrt(N); +A = A(1: 256, :); + +Target_index = 321; +SNR = 10; +lambda = 0.0005; + +%% +[ x_mf, x_cs ] = yAxn_recovery( A, SNR, Target_index, lambda ); + +figure(1111) +subplot(211) +plot(abs(x_mf)) +subplot(212) +plot(abs(x_cs)) + + +%% +rep_time = 500; +[ res_MF, res_CS, P_fa, H1_MF_cnt, H0_MF_cnt, H1_CS_cnt, H0_CS_cnt ] = yAxn_MC( A, SNR, Target_index, lambda, rep_time ); + +P_fa_MF = res_MF(:,1); +P_d_MF = res_MF(:,2); +P_fa_CS = res_CS(:,1); +P_d_CS = res_CS(:,2); + +figure; +loglog(P_fa,P_fa_MF, 'linewidth', 2); +hold on; +loglog(P_fa,P_fa_CS, 'linewidth', 2); +legend('MF','CS') +xlabel('P_{fa}'); +ylabel('Actual P_{fa}'); + +figure; +semilogx(P_fa,P_d_MF, 'linewidth', 2); +hold on; +semilogx(P_fa,P_d_CS, 'linewidth', 2); +legend('MF','CS') +xlabel('P_{fa} set'); +ylabel('P_{d}'); + +figure; +semilogx(P_fa_MF,P_d_MF, 'linewidth', 2); +hold on; +grid on; +semilogx(P_fa_CS,P_d_CS, 'linewidth', 2); +xlim([1e-4,1]) +legend('MF','CS') +xlabel('Actual P_{fa}'); +ylabel('P_{d}'); +title('ROC'); + + +figure +subplot(211) +histogram(real(H0_MF_cnt(5,:))) +hold on; +histogram(real(H1_MF_cnt)) +subplot(212) +histogram(real(H0_CS_cnt(5,:))) +hold on; +histogram(real(H1_CS_cnt)) diff --git a/0 - example/wzk - PSK_ROC/test_signal_t.m b/0 - example/wzk - PSK_ROC/test_signal_t.m new file mode 100644 index 0000000..cd67a33 --- /dev/null +++ b/0 - example/wzk - PSK_ROC/test_signal_t.m @@ -0,0 +1,57 @@ +clc +clear +close all + +rng(6) + +PRF = 5000; % 脉冲重复频率 +Tr = 1 / PRF; % 脉冲重复间隔 +Tp = 2e-5; +P = 10; +B = P / Tp; +fs = 10 * B; +Ts = 1/fs; +fc = 1.25e9; + +t = 0:1/fs:1-1/fs; + +tsin = Tp / P * 2; +fsin = 1 / tsin; + +N = round(Tr * fs); +N_high = round(Tp * fs / P); + +signal_t = zeros(1, N); +signal_td = zeros(1, N); + +sign_p = sign(randn(1, P)); + +for p = 1: P + for i = 1: N_high + tp = (i - 1) * (1 / fs); + signal_t(1, (p-1)*N_high + i) = sign_p(p) * exp(1j * 2 * pi * fsin * tp); + tp2 = (i - 0.5) * (1 / fs); + signal_td(1, (p-1)*N_high + i) = sign_p(p) * exp(1j * 2 * pi * fsin * tp2); + end +end + +A = zeros(N, 2 * N); +A(:, 1) = transpose(signal_t); +A(:, 2) = transpose([0, signal_td(1: end-1)]); +for k = 3: 2 * N + A(:, k) = [0; A(1: end - 1, k - 2)]; +end + + + +figure(1) +subplot(211) +plot(real(A(1:N_high*P+20, 1))) +hold on; +plot(real(A(1:N_high*P+20, 2))) +plot(real(A(1:N_high*P+20, 3))) +subplot(212) +plot(imag(A(1:N_high*P+20, 1))) +hold on; +plot(imag(A(1:N_high*P+20, 2))) +plot(imag(A(1:N_high*P+20, 3))) diff --git a/0 - example/wzk - PSK_ROC/yAxn_MC.m b/0 - example/wzk - PSK_ROC/yAxn_MC.m new file mode 100644 index 0000000..049d163 --- /dev/null +++ b/0 - example/wzk - PSK_ROC/yAxn_MC.m @@ -0,0 +1,106 @@ +function [ res_MF, res_CS, P_fa, H1_MF_cnt, H0_MF_cnt, H1_CS_cnt, H0_CS_cnt ] = yAxn_MC( A, SNR, Target_index, lambda, rep_time ) + + % matrix parameters + [M, N] = size(A); + mul = A(:,1)'*A(:,1); + + % H0 sample + L = round(Target_index + N * 0.1); + R = round(N * 0.9); + Lambda_C = L: R; + + % parameters + sigma_n = 0.1; + Target_amplitude = sqrt(10^(SNR/10) * sigma_n^2 / mul); + + + %% CS setting + % CS-parameters +% lambda = 0.0005; + alpha = 1/4; + delta = 1e-8*alpha; + + % CS-normalization + J1 = A*A'; + lambda_J=eig(J1); + A_norm = A / sqrt(lambda_J(end)); + Hp = A_norm'*A_norm; + [~, D] = eig(Hp); + d = diag(D); + sigma_n_norm = sigma_n / sqrt(lambda_J(end)); + + + %% generate x + x = zeros(N, 1); + x(Target_index) = Target_amplitude; + + + %% MC parameters + P_fa = [1e-6, 1e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1]; + len_P_fa = length(P_fa); + + P_fa_MF_cnt = zeros(len_P_fa, rep_time); + P_d_MF_cnt = zeros(len_P_fa, rep_time); + P_fa_CS_cnt = zeros(len_P_fa, rep_time); + P_d_CS_cnt = zeros(len_P_fa, rep_time); + + H1_index = zeros(N, 1); + H1_index(Target_index) = 1; + H0_index = ones(size(x)); + H0_index(Target_index) = 0; + H1_MF_cnt = zeros(1, rep_time); + H0_MF_cnt = zeros(N - 1, rep_time); + H1_CS_cnt = zeros(1, rep_time); + H0_CS_cnt = zeros(N - 1, rep_time); + + %% MC + parfor rep = 1: rep_time +% for rep = 1: rep_time + %% generate y + noise = random('Normal', 0, sigma_n/sqrt(2), M, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), M, 1); + y = A*x + noise; + + %% recover + % MF + x_mf = A' * y ./ mul; + % sigma_MF = sqrt(var(x_mf(Lambda_C))); + sigma_MF = sigma_n / sqrt(mul); + x_mf_norm = x_mf ./ sigma_MF; + stat_MF = abs(x_mf ./ sigma_MF).^2; + + % CS + y_norm = y / sqrt(lambda_J(end)); + x_FISTA = FISTA_v1(y_norm, A_norm, lambda, delta, Hp); + [x_d_cal_f, hat_Q1_cal_f, sigma_d_cal_f] = cal_debiased_LASSO_v1(x_FISTA, A_norm, y_norm, lambda, sigma_n_norm, d); + % sigma_CS = sqrt(var(x_d_cal_f(Lambda_C))); + sigma_CS = sigma_d_cal_f; + x_cs_norm = x_d_cal_f ./ sigma_CS; + stat_CS = abs(x_d_cal_f ./ sigma_CS).^2; + + %% detect + H1_MF_cnt(rep) = x_mf(Target_index); + H0_MF_cnt(:,rep) = [x_mf(1:Target_index-1);x_mf(Target_index+1:end)]; + H1_CS_cnt(rep) = x_d_cal_f(Target_index); + H0_CS_cnt(:,rep) = [x_d_cal_f(1:Target_index-1);x_d_cal_f(Target_index+1:end)]; + + kd = chi2inv(1 - P_fa, 2) / 2; + for cnt_h_th = 1: len_P_fa + P_fa_MF_cnt(cnt_h_th, rep) = sum(stat_MF(H0_index > 0) > kd(cnt_h_th)) / sum(H0_index); + P_d_MF_cnt(cnt_h_th, rep) = sum(stat_MF(H1_index > 0) > kd(cnt_h_th)) / sum(H1_index); + P_fa_CS_cnt(cnt_h_th, rep) = sum(stat_CS(H0_index > 0) > kd(cnt_h_th)) / sum(H0_index); + P_d_CS_cnt(cnt_h_th, rep) = sum(stat_CS(H1_index > 0) > kd(cnt_h_th)) / sum(H1_index); + end + + fprintf('%d\n', rep); + end + + P_fa_MF = mean(P_fa_MF_cnt, 2); + P_d_MF = mean(P_d_MF_cnt, 2); + P_fa_CS = mean(P_fa_CS_cnt, 2); + P_d_CS = mean(P_d_CS_cnt, 2); + + res_MF = [P_fa_MF, P_d_MF]; + res_CS = [P_fa_CS, P_d_CS]; + +end + diff --git a/0 - example/wzk - PSK_ROC/yAxn_recovery.m b/0 - example/wzk - PSK_ROC/yAxn_recovery.m new file mode 100644 index 0000000..508c347 --- /dev/null +++ b/0 - example/wzk - PSK_ROC/yAxn_recovery.m @@ -0,0 +1,57 @@ +function [ x_mf_norm, x_cs_norm ] = yAxn_recovery( A, SNR, Target_index, lambda ) + + % matrix parameters + [M, N] = size(A); + mul = A(:,1)'*A(:,1); + + % H0 sample + L = round(Target_index + N * 0.1); + R = round(N * 0.9); + Lambda_C = L: R; + + % parameters + sigma_n = 0.1; + Target_amplitude = sqrt(10^(SNR/10) * sigma_n^2 / mul); + + + %% CS setting + % CS-parameters +% lambda = 0.0005; + alpha = 1/4; + delta = 1e-8*alpha; + + % CS-normalization + J1 = A*A'; + lambda_J=eig(J1); + A_norm = A / sqrt(lambda_J(end)); + Hp = A_norm'*A_norm; + [~, D] = eig(Hp); + d = diag(D); + sigma_n_norm = sigma_n / sqrt(lambda_J(end)); + + + %% generate x + x = zeros(N, 1); + x(Target_index) = Target_amplitude; + + + %% generate y + noise = random('Normal', 0, sigma_n/sqrt(2), M, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), M, 1); + y = A*x + noise; + + + %% recover + % MF + x_mf = A' * y ./ mul; + sigma_MF = sqrt(var(x_mf(Lambda_C))); + x_mf_norm = x_mf ./ sigma_MF; + + % CS + y_norm = y / sqrt(lambda_J(end)); + x_FISTA = FISTA_v1(y_norm, A_norm, lambda, delta, Hp); + [x_d_cal_f, hat_Q1_cal_f, sigma_d_cal_f] = cal_debiased_LASSO_v1(x_FISTA, A_norm, y_norm, lambda, sigma_n_norm, d); + sigma_CS = sqrt(var(x_d_cal_f(Lambda_C))); + x_cs_norm = x_d_cal_f ./ sigma_CS; + +end + diff --git a/0 - example/wzk - RDA/angle_estimation.m b/0 - example/wzk - RDA/angle_estimation.m new file mode 100644 index 0000000..7f1a0e3 --- /dev/null +++ b/0 - example/wzk - RDA/angle_estimation.m @@ -0,0 +1,60 @@ +function [ target_list_RDA ] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda, acc) +% This function estimates all targets angle +% +% Usage: +% [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda, acc) +% +% Inputs: +% target_list_RDA_temp: Temporary list of targets before angle estimation +% echo_mtx: Original echo data matrix +% A: Chirp measurement matrix +% d: Antenna spacing +% lambda: Wave length +% acc: Accuracy for angle estimation +% +% Outputs: +% target_list_RDA: Target list include range, doppler and angle information + + + % parameters + if nargin < 6 + acc = 181; + end + + % calculate all targets angle + target_list_RDA = zeros(size(target_list_RDA_temp)); + for i = 1: size(target_list_RDA_temp, 1) + t_rda = target_list_RDA_temp(i, :); + target_list_RDA(i, 1:2) = t_rda(1:2); + target_list_RDA(i, 3) = cal_angle(t_rda(1), t_rda(2), echo_mtx, A, d, lambda, acc); + end + +end + + +%% sub-functions +% calculate target actual angle with range and doppler index +function [ angle ] = cal_angle(r_idx, d_idx, echo_mtx, A, d, lambda, acc) + + % range doppler processing + echo_r = zeros([size(echo_mtx, 1), 1, size(echo_mtx, 3)]); + for i = 1 : size(echo_mtx, 3) + echo_r(:, 1, i) = sum(bsxfun(@times, echo_mtx(:, :, i), A(:, r_idx)'), 2); + end + echo_r = squeeze(echo_r); + echo_r_ex = [echo_r, zeros(size(echo_r))]; + echo_rd = zeros(size(echo_r_ex, 1), 1); + for i = 1 : size(echo_r_ex, 1) + temp_v = fftshift(fft(echo_r_ex(i,:))); + echo_rd(i) = temp_v(d_idx); + end + + % CBF + THETA = linspace(-90, 90, round(acc)); + fai = exp((0: size(echo_r_ex, 1) - 1)' *... + -1j * 2 * pi / lambda * d * sin(THETA / 180 * pi)); + angle_mf = abs(fai'* echo_rd); + [v, idx] = max(angle_mf); + angle = THETA(idx); + +end diff --git a/0 - example/wzk - RDA/doppler_process_CS.m b/0 - example/wzk - RDA/doppler_process_CS.m new file mode 100644 index 0000000..33c11a3 --- /dev/null +++ b/0 - example/wzk - RDA/doppler_process_CS.m @@ -0,0 +1,200 @@ +function [ echo_rd_mtx, stat_RD, sigma_n_o ] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta) +% Doppler processing for echo data using compressed sensing +% +% Usage: +% [echo_rd_mtx, stat_RD, sigma_n_o] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta) +% +% Inputs: +% echo_r_mtx: Range processing echo data +% sigma_n: Input noise standard deviation +% lambda: LASSO weight +% gamma: Compressed ratio(0.5) +% delta: Convergence normalized difference(1e-6) +% +% Outputs: +% echo_rd_mtx: Range-doppler processing echo data +% stat_RD: Range-doppler statistics +% sigma_n_o: Output noise standard deviation + + + % parameters + if nargin < 5 + delta = 1e-6; + end + if nargin < 4 + gamma = 0.5; + end + [lenA, lenR, M] = size(echo_r_mtx); + N = round(M / gamma); + iter_max_VAMP = 1000; + lambda_v = zeros(N, 1) + lambda; + + % generate mtx + F_ori = dftmtx(N); + F = F_ori(1:M,:); + F_inv = conj(F) / N; + + % normalization + A = (sqrt(N) * eye(M)) * F_inv; + echo_r_mtx = sqrt(N) .* echo_r_mtx; + sigma_n = sqrt(N) * sigma_n; + + % doppler processing + echo_rd_mtx = zeros(lenA, lenR, N); + stat_RD = zeros(lenA, lenR, N); + sigma_n_o_cnt = zeros(lenA, lenR); + for numA = 1: lenA + parfor numR = 1: lenR + sample = squeeze(echo_r_mtx(numA, numR, :)); + y = sample; + x_LASSO = cVAMPro(y, A, lambda_v, delta, iter_max_VAMP); + [x_d_CROD, sigma_CROD] = CROD(y, A, x_LASSO, lambda, sigma_n); + sigma_n_o_cnt(numA, numR) = abs(sigma_CROD); + stat_RD(numA, numR, :) = abs(fftshift(x_d_CROD) / sigma_CROD).^2; + echo_rd_mtx(numA, numR, :) = fftshift(x_d_CROD); + end + end +% sigma_n_o = mean(sigma_n_o_cnt, 2); + sigma_n_o = mean(mean(sigma_n_o_cnt)); + +end + + +%% sub-functions +% algorithm for LASSO +% y: measurements +% A: measurement matrix +% lambda: LASSO weight +% tau: convergence normalized difference +% Kit: maximum number of iterations +% LASSO estimator +function x_hat_wl = cVAMPro(y, A, lambda, tau, Kit) + + % Initialization + [M, N] = size(A); + gamma = M / N; + k = 0; + p = ctranspose(A) * y; + h_1 = p; + Q_1 = gamma; + tau_d = 1; + + % Iteration + while ((k < Kit) && (tau_d > tau)) + % Factorized Part + x_1 = ST(h_1, lambda, Q_1); + chi_1 = F1(x_1, lambda, Q_1); + % Message Passing + h_2 = x_1 / chi_1 - h_1; + Q_2 = 1 / chi_1 - Q_1; + % Gaussian Part + t1 = (p + h_2) / Q_2; + t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2); + x_2 = t1 - t2; + chi_2 = gamma / (Q_2 + 1) + (1 - gamma) / Q_2; + % Message Passing + h_1_next = x_2 ./ chi_2 - h_2; + Q_1_next = 1 / chi_2 - Q_2; + tau_d = norm(h_1_next - h_1, Inf) / norm(h_1_next, Inf); + k = k + 1; + % output + x_hat_wl = x_1; + % next + h_1 = h_1_next; + Q_1 = Q_1_next; + end + +end + + +% soft threshold function +% x: processing object +% thd: threshold +% y: result +function x = ST(h_1, lambda, Q_1) + [N, M] = size(h_1); + x = zeros(N, M); + + for i = 1:N + sign = h_1(i) ./ abs(h_1(i)); + diff = abs(h_1(i)) - lambda(i); + x(i) = sign .* (diff ./ Q_1) .* SF(diff); + end + +end + + +% Heaviside's step function +function v = SF(a) + + if a > 0 + v = 1; + elseif a == 0 + v = 0; % at zero points + else + v = 0; + end + +end + + +% Calculation of chi_1 +function chi_1 = F1(x_1, lambda, Q_1) + + [N, M] = size(x_1); + count = 0; + + for i = 1:N + temp = Q_1 * abs(x_1(i)) + lambda(i); + count = count + (2 - lambda(i) / temp) * SF(abs(x_1(i))); +% count = count + (2-lambda(i)/temp) * (abs(x_1(i)) > 1e-4); + end + + chi_1 = count / (2 * N * Q_1); + +end + +% calculate debiased LASSO estimator +% y: measurements +% A: measurement matrix +% x_LASSO: LASSO estimator +% lambda: LASSO weight +% sigma_n: input noise standard deviation +% x_d_CROD: debiased LASSO estimator +% sigma_CROD: equivalent noise standard deviation estimator +function [ x_d_CROD, sigma_CROD ] = CROD(y, A, x_LASSO, lambda, sigma_n) + + [m, n] = size(A); + gamma = m / n; + + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat; + + RSS = sum(abs(y - A * x_LASSO).^2)/m; + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + +end diff --git a/0 - example/wzk - RDA/doppler_process_MF.m b/0 - example/wzk - RDA/doppler_process_MF.m new file mode 100644 index 0000000..3d7b0ce --- /dev/null +++ b/0 - example/wzk - RDA/doppler_process_MF.m @@ -0,0 +1,45 @@ +function [ echo_rd_mtx, stat_RD, sigma_n_o ] = doppler_process_MF(echo_r_mtx, sigma_n, gamma) +% Doppler processing for echo data using matching filter +% +% Usage: +% [echo_rd_mtx, stat_RD, sigma_n_o] = doppler_process_MF(echo_r_mtx, sigma_n, gamma) +% +% Inputs: +% echo_r_mtx: Range processing echo data +% sigma_n: Input noise standard deviation +% gamma: Compressed ratio(0.5) +% +% Outputs: +% echo_rd_mtx: Range-doppler processing echo data +% stat_RD: Range-doppler statistics +% sigma_n_o: Output noise standard deviation + + + % parameters + if nargin < 3 + gamma = 0.5; + end + [lenA, lenR, M] = size(echo_r_mtx); + N = round(M / gamma); + + % generate mtx + F_ori = dftmtx(N); + F = F_ori(1:M,:); + multiple_d = F(:,1)' * F(:,1); + + % doppler matched filtering + echo_rd_mtx = zeros(lenA, lenR, N); + for numA = 1: lenA + for numR = 1: lenR + sample = squeeze(echo_r_mtx(numA, numR, :)); + dpl_temp = transpose(F) * sample; + dpl_norm = fftshift(dpl_temp ./ multiple_d); + echo_rd_mtx(numA, numR, :) = dpl_norm; + end + end + + % calculate output noise + sigma_n_o = sqrt(sigma_n^2 / multiple_d); + stat_RD = abs(echo_rd_mtx ./ sigma_n_o).^2; + +end diff --git a/0 - example/wzk - RDA/echo_processing_CS.m b/0 - example/wzk - RDA/echo_processing_CS.m new file mode 100644 index 0000000..30ea5f7 --- /dev/null +++ b/0 - example/wzk - RDA/echo_processing_CS.m @@ -0,0 +1,56 @@ +function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] = ... + echo_processing_CS( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d ) +% This function processes echo data using compressed sensing method. +% +% Usage: +% [echo_rd_mtx, target_list_RDA, sigma_n_rd] = echo_processing_MF(echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d) +% +% Inputs: +% echo_mtx: Matrix containing the echo data +% PRF: Pulse Repetition Frequency +% fs: Sampling frequency +% fc: Carrier frequency +% B: Bandwidth +% D: Duty ratio +% d: Antenna spacing +% sigma_n: Input noise level of the echo data +% P_fa: False alarm probability threshold +% lambda_r: LASSO weight for range processing +% lambda_d: LASSO weight for doppler processing +% +% Outputs: +% echo_rd_mtx: Echo matrix after range-Doppler processing +% target_list_RDA: List of detected targets +% sigma_n_rd: Estimated noise level after range-Doppler processing + + + % parameters + c = 3e8; + lambda = c / fc; + [num_antenna, N, num_pluse] = size(echo_mtx); + + % data processing + [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse); + + [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda); + fprintf(' Spatial_FFT done.\n'); + + [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + fprintf(' Generate_chirp_mtx done.\n'); + + [echo_r_mtx, sigma_n_o] = range_process_CS(echo_sFFT_mtx, A, sigma_n, lambda_r); + fprintf(' Range_processing done.\n'); + + [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_CS(echo_r_mtx, sigma_n_o, lambda_d); + fprintf(' Doppler_processing done.\n'); + + [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa); + fprintf(' Target_detection done.\n'); + + [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda); + fprintf(' Angle_estimation done.\n'); + target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1)); + target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2)); + +end + diff --git a/0 - example/wzk - RDA/echo_processing_CS_v0.m b/0 - example/wzk - RDA/echo_processing_CS_v0.m new file mode 100644 index 0000000..08c9710 --- /dev/null +++ b/0 - example/wzk - RDA/echo_processing_CS_v0.m @@ -0,0 +1,28 @@ +function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] = ... + echo_processing_CS_v0( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d ) + + % parameters + c = 3e8; + lambda = c / fc; + [num_antenna, N, num_pluse] = size(echo_mtx); + + + % data processing + [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse); + + [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda); + + [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + + [echo_r_mtx, sigma_n_o] = range_process_CS(echo_sFFT_mtx, A, sigma_n, lambda_r); + + [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_CS(echo_r_mtx, sigma_n_o, lambda_d); %差个方差 + + [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa); + target_list_RDA = target_list_RDA_temp; + target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1)); + target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2)); + target_list_RDA(:, 3) = angle_v(target_list_RDA(:, 3)); + +end + diff --git a/0 - example/wzk - RDA/echo_processing_MF.m b/0 - example/wzk - RDA/echo_processing_MF.m new file mode 100644 index 0000000..3d69dec --- /dev/null +++ b/0 - example/wzk - RDA/echo_processing_MF.m @@ -0,0 +1,54 @@ +function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] =... + echo_processing_MF( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa ) +% This function processes echo data using matching filter method. +% +% Usage: +% [echo_rd_mtx, target_list_RDA, sigma_n_rd] = echo_processing_MF(echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa) +% Inputs: +% echo_mtx: Matrix containing the echo data +% PRF: Pulse Repetition Frequency +% fs: Sampling frequency +% fc: Carrier frequency +% B: Bandwidth +% D: Duty ratio +% d: Antenna spacing +% sigma_n: Input noise level of the echo data +% P_fa: False alarm probability threshold +% +% Outputs: +% echo_rd_mtx: Echo matrix after range-Doppler processing +% target_list_RDA: List of detected targets +% sigma_n_rd: Estimated noise level after range-Doppler processing + + + % parameters + c = 3e8; + lambda = c / fc; + [num_antenna, N, num_pluse] = size(echo_mtx); + + + % data processing + [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse); + + [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda); + fprintf(' Spatial_FFT done.\n'); + + [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + fprintf(' Generate_chirp_mtx done.\n'); + + [echo_r_mtx, sigma_n_o] = range_process_MF(echo_sFFT_mtx, A, sigma_n); + fprintf(' Range_processing done.\n'); + + [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_MF(echo_r_mtx, sigma_n_o); + fprintf(' Doppler_processing done.\n'); + + [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa); + fprintf(' Target_detection done.\n'); + + [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda); + fprintf(' Angle_estimation done.\n'); + target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1)); + target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2)); + +end + diff --git a/0 - example/wzk - RDA/echo_processing_MF_v0.m b/0 - example/wzk - RDA/echo_processing_MF_v0.m new file mode 100644 index 0000000..8b0085a --- /dev/null +++ b/0 - example/wzk - RDA/echo_processing_MF_v0.m @@ -0,0 +1,28 @@ +function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] =... + echo_processing_MF_v0( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa ) + + % parameters + c = 3e8; + lambda = c / fc; + [num_antenna, N, num_pluse] = size(echo_mtx); + + + % data processing + [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse); + + [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda); + + [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + + [echo_r_mtx, sigma_n_o] = range_process_MF(echo_sFFT_mtx, A, sigma_n); + + [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_MF(echo_r_mtx, sigma_n_o); %差个方差 + + [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa); + target_list_RDA = target_list_RDA_temp; + target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1)); + target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2)); + target_list_RDA(:, 3) = angle_v(target_list_RDA(:, 3)); + +end + diff --git a/0 - example/wzk - RDA/generate_chirp_mtx.m b/0 - example/wzk - RDA/generate_chirp_mtx.m new file mode 100644 index 0000000..993bd66 --- /dev/null +++ b/0 - example/wzk - RDA/generate_chirp_mtx.m @@ -0,0 +1,81 @@ +function [ A, signal_t ] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid) +% Generates a chirp measurement matrix +% +% Usage: +% [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid) +% +% Inputs: +% PRF: Pulse Repetition Frequency +% B: Bandwidth +% fs: Sampling frequency +% D: Duty ratio +% gamma: Compression ratio +% sign_mid: if sign_mid is 0 means the initial frequency is 0, +% if sign_mid is 1 means the initial frequency is -B/2, +% +% Outputs: +% A: Chirp measurement matrix +% signal_t: Transmitting signal + + + % parameters + if nargin < 5 + gamma = 0.5; + end + if nargin < 6 + sign_mid = 0; + end + Tr = 1 / PRF; + Tp = Tr * D; + K = B / Tp; + + % generate_signal + N = Tr * fs; + N_high = Tp * fs; + N_mtx = round(N / gamma); + + signal_t = zeros(1, N); + for i = 1: N_high + tp = i * (1 / fs) - sign_mid * N_high / fs / 2; + signal_t(1, i) = exp(1j * 2 * pi * 0.5 * K * tp .^ 2); + end + + signal_t_2fs = zeros(1, N_mtx); + for i = 1: round(N_high / gamma) + tp = (i+1) * (1 / fs / 2) - sign_mid * N_high / fs / 2; + signal_t_2fs(1, i) = exp(1j * pi * K * tp .^ 2); + end + + % generate chirp matrix + A = generate_matrix_by_signal2fs(transpose(signal_t_2fs)); + +end + + +%% sub-functions +% generate chirp matrix when gamma=0.5 +function [ mtx ] = generate_matrix_by_signal2fs(signal) + + mtx = []; + l = round(length(signal) / 2); + temp1 = signal(1:2:end); + temp2 = circshift(signal(2:2:end), 1); + if length(temp2) ~= length(temp1) + temp2 = [temp2; 0]; + end + + for i = 1:l + t1 = circshift(temp1, i-1); + t2 = circshift(temp2, i-1); + if i - 1 > 0 + t1(1:i - 1,1) = 0; + end + if i - 1 > 0 + t2(1:i - 1,1) = 0; + end + mtx = [mtx,t1,t2]; + end + +end + + diff --git a/0 - example/wzk - RDA/get_range_speed_val.m b/0 - example/wzk - RDA/get_range_speed_val.m new file mode 100644 index 0000000..ea7de79 --- /dev/null +++ b/0 - example/wzk - RDA/get_range_speed_val.m @@ -0,0 +1,35 @@ +function [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, numP, gamma_r, gamma_d) +% Calculates range and speed values for node index +% +% Usage: +% [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, numP, gamma_r, gamma_d) +% +% Inputs: +% PRF: Pulse Repetition frequency +% fs: Sampling frequency +% fc: Carrier frequency +% numP: Number of pulses +% gamma_r: Range compression ratio(0.5) +% gamma_d: Doppler compression ratio(0.5) +% +% Outputs: +% distance_v: Real distance value +% speed_v: Real speed value + + + % parameters + if nargin < 6 + gamma_d = 0.5; + end + if nargin < 5 + gamma_r = 0.5; + end + c = 3e8; + Tr = 1 / PRF; + lambda = c / fc; + + % calculate + distance_v = 0: (c / fs / 2 * gamma_r) : (Tr * c / 2 - c / fs / 2 * gamma_r); + speed_v = -(-PRF / 2: PRF / numP * gamma_d: PRF / 2 - PRF / numP * gamma_d) * lambda / 2 ; + +end \ No newline at end of file diff --git a/0 - example/wzk - RDA/main_CS.m b/0 - example/wzk - RDA/main_CS.m new file mode 100644 index 0000000..fb954ee --- /dev/null +++ b/0 - example/wzk - RDA/main_CS.m @@ -0,0 +1,81 @@ +clc +clear +close all + +%% load echo data +filename = 'Raw_Echo_5dB'; +load(['./data/', filename, '.mat']); + + +%% parameters +% ladar +PRF = 5000; +B = 5e6; +D = 0.1; +Tp = 2e-5; +fs = 5e6; +fc = 1.25e9; + +% antenna +num_antenna = 18; +d = 0.12; + +% echo +num_pulse = 64; +sigma_n = 0.1; + +% detect +P_fa = 1e-6; + +% algorithm +gamma_r = 0.5; +gamma_d = 0.5; +lambda_r = 0.005; +lambda_d = 0.15; + +% accurate angle +sign_AA = 1; + + +%% data processing +fprintf('Using CS\n'); +fprintf(['File: ', filename, '\n']); + +tic; +if sign_AA == 0 + [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ... + echo_processing_CS_v0( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d ); +else + [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ... + echo_processing_CS( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d ); +end +toc; + + +%% plot +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + +figure(1) +scatter3(target_list_RDA(:,1), target_list_RDA(:,2),... + target_list_RDA(:,3), 'filled', 'o') +xlabel('Range(m)'); +ylabel('Speed(m/s)'); +zlabel('Angle(°)'); +zlim([-90, 90]); +ylim([80, 150]); +xlim([17000, 19000]); +set(gca, 'FontSize', Fontsize); +title('node') +set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]); +set(gca,'fontsize',18,'fontname','Times'); + + +%% save +save(['./output/', filename, '_CS_Pfa', num2str(P_fa), '.mat'],... + 'target_list_RDA', 'echo_rd_mtx', 'sigma_n_out', 'P_fa'); +fprintf(['[', filename, '_CS_Pfa', num2str(P_fa), '.mat]', ' saved.\n\n']); + diff --git a/0 - example/wzk - RDA/main_MF.m b/0 - example/wzk - RDA/main_MF.m new file mode 100644 index 0000000..fafd12a --- /dev/null +++ b/0 - example/wzk - RDA/main_MF.m @@ -0,0 +1,79 @@ +clc +clear +close all + +%% load echo data +filename = 'Raw_Echo_5dB'; +load(['./data/', filename, '.mat']); + + +%% parameters +% ladar +PRF = 5000; +B = 5e6; +D = 0.1; +Tp = 2e-5; +fs = 5e6; +fc = 1.25e9; + +% antenna +num_antenna = 18; +d = 0.12; + +% echo +num_pulse = 64; +sigma_n = 0.1; + +% detect +P_fa = 1e-6; + +% algorithm +gamma_r = 0.5; +gamma_d = 0.5; + +% accurate angle +sign_AA = 1; + + +%% data processing +fprintf('Using MF\n'); +fprintf(['File: ', filename, '\n']); + +tic; +if sign_AA == 0 + [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ... + echo_processing_MF_v0( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa ); +else + [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ... + echo_processing_MF( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa ); +end +toc; + + +%% plot +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + +figure(1) +scatter3(target_list_RDA(:,1), target_list_RDA(:,2),... + target_list_RDA(:,3), 'filled', 'o') +xlabel('Range(m)'); +ylabel('Speed(m/s)'); +zlabel('Angle(°)'); +zlim([-90, 90]); +ylim([80, 150]); +xlim([17000, 19000]); +set(gca, 'FontSize', Fontsize); +title('node') +set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]); +set(gca,'fontsize',18,'fontname','Times'); + + +%% save +save(['./output/', filename, '_MF_Pfa', num2str(P_fa), '.mat'],... + 'target_list_RDA', 'echo_rd_mtx', 'sigma_n_out', 'P_fa'); +fprintf(['[', filename, '_MF_Pfa', num2str(P_fa), '.mat]', ' saved.\n\n']); + diff --git a/0 - example/wzk - RDA/plot_result.m b/0 - example/wzk - RDA/plot_result.m new file mode 100644 index 0000000..6d810e8 --- /dev/null +++ b/0 - example/wzk - RDA/plot_result.m @@ -0,0 +1,38 @@ +clc +clear +close all + + +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + + +Pfa_set = '1e-06'; +method = 'CS'; +SNR = 5; + +filename = ['Raw_Echo_', num2str(SNR), 'dB_', method, '_Pfa', Pfa_set]; + +load(['./output/', filename, '.mat']); + +%% +figure(1) +scatter3(target_list_RDA(:,1), target_list_RDA(:,2),... + target_list_RDA(:,3), 'filled', 'o') +xlabel('Range(m)'); +ylabel('Speed(m/s)'); +zlabel('Angle(°)'); +zlim([-90, 90]); +% ylim([80, 150]); +% xlim([17000, 19000]); +set(gca, 'FontSize', Fontsize); +title('node') +set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]); +set(gca,'fontsize',18,'fontname','Times'); + + + + diff --git a/0 - example/wzk - RDA/range_process_CS.m b/0 - example/wzk - RDA/range_process_CS.m new file mode 100644 index 0000000..fa565dd --- /dev/null +++ b/0 - example/wzk - RDA/range_process_CS.m @@ -0,0 +1,155 @@ +function [ echo_r_mtx, sigma_n_o ] = range_process_CS(echo_mtx, A, sigma_n, lambda, delta) +% Range processing for echo data using compressed sensing +% +% Usage: +% [echo_r_mtx, sigma_n_o] = range_process_CS(echo_mtx, A, sigma_n, lambda, delta) +% +% Inputs: +% echo_mtx: Original echo data +% A: Chirp measurement matrix +% sigma_n: Input noise standard deviation +% lambda: LASSO weight +% delta: Convergence normalized difference(2e-7) +% +% Outputs: +% echo_r_mtx: Range processing echo data +% sigma_n_o: Output noise standard deviation + + + % parameters + if nargin < 5 + delta = 2e-7; + end + [M, N] = size(A); + [lenA, M, lenP] = size(echo_mtx); + + % normalization + J1 = A*A'; + lambda_J=eig(J1); + A = A / sqrt(lambda_J(end)); + echo_mtx = echo_mtx ./ sqrt(lambda_J(end)); + sigma_n = sigma_n / sqrt(lambda_J(end)); + + % compressed sensing + echo_r_mtx = zeros(lenA, N, lenP); + sigma_n_o_cnt = zeros(lenA, lenP); + for numA = 1: lenA + parfor numP = 1: lenP + sr = echo_mtx(numA, :, numP); + y = transpose(sr); + % LASSO + x_FISTA = FISTA(y, A, lambda, delta); + % debiased LASSO + [x_d, sigma_w] = cal_debiased_LASSO(x_FISTA, A, y, lambda, sigma_n); + sigma_n_o_cnt(numA, numP) = abs(sigma_w); + echo_r_mtx(numA, :, numP) = x_d; + end + end + + % calculate output noise + sigma_n_o = mean(mean(sigma_n_o_cnt)); + +end + + +%% sub-functions +% algorithm for LASSO +% y: measurements +% A: measurement matrix +% lambda: LASSO weight +% delta: convergence normalized difference +% z: LASSO estimator +function [z] = FISTA(y, A, lambda, delta) + + x_pre = A'*y; + t = 1; + z = x_pre; + z_pre = z; + t_pre = t; + N = size(A, 2); + diff = 1; + E = eig(A'*A); + L = E(end); + temp1 = A'*y/L; + temp2 = eye(N) - A'*A/L; + k = 0; + + while((diff > delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = mean(abs(z_pre - z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; + end + +end + + +% soft threshold function +% x: processing object +% thd: threshold +% y: result +function y = sft_thd(x, thd) + + if isequal(size(x), size(thd)) + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd); + else + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd); + end + +end + + +% calculate debiased LASSO estimator +% x: LASSO estimator +% A: measurement matrix +% y: measurements +% lambda: LASSO weight +% sigma_n: input noise standard deviation +% x_d: debiased LASSO estimator +% sigma_d: equivalent noise standard deviation +function [x_d, sigma_d] = cal_debiased_LASSO(x, A, y, lambda, sigma) + + [M, N] = size(A); + gamma = M/N; + hat_Q1 = gamma; + [~, D] = eig(A'*A); + d = diag(D); + + diff = 1; + T = 1000; + t = 0; + + while (t < T) && (diff > 1e-6) + Q1_pre = hat_Q1; + rho = mean((2 - lambda./(hat_Q1*abs(x) + lambda)).*(abs(x) > 1e-4))/2; + hat_Q1 = rho/mean(1./(d + (1-rho)*hat_Q1/rho)); + diff = abs(Q1_pre - hat_Q1); + t = t+1; + end + + x_d = x + 1/hat_Q1*A'*(y - A*x); + + chi = rho/hat_Q1; + hat_Q2 = 1/chi - hat_Q1; + + t = -hat_Q2; + t_prime = -1/mean((1./(d+hat_Q2)).^2); + G_prime = t + 1/chi; + G_wprime = t_prime + 1/chi/chi; + RSS = sum(abs(y - A*x).^2)/M; + hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +... + (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2; + sigma_d = sqrt(2*hat_chi)/hat_Q1; + +end diff --git a/0 - example/wzk - RDA/range_process_MF.m b/0 - example/wzk - RDA/range_process_MF.m new file mode 100644 index 0000000..3801e43 --- /dev/null +++ b/0 - example/wzk - RDA/range_process_MF.m @@ -0,0 +1,36 @@ +function [ echo_r_mtx, sigma_n_o ] = range_process_MF(echo_mtx, A, sigma_n) +% Range processing for echo data using matching filter +% +% Usage: +% [echo_r_mtx, sigma_n_o] = range_process_MF(echo_mtx, A, sigma_n) +% +% Inputs: +% echo_mtx: Original echo data +% A: Chirp measurement matrix +% sigma_n: Input noise standard deviation +% +% Outputs: +% echo_r_mtx: Range processing echo data +% sigma_n_o: Output noise standard deviation + + + % parameters + multiple_r = A(:, 1)' * A(:, 1); + [M, N] = size(A); + [lenA, M, lenP] = size(echo_mtx); + + % matched filtering + echo_r_mtx = zeros(lenA, N, lenP); + for numA = 1: lenA + for numP = 1: lenP + sr = echo_mtx(numA, :, numP); + mf_res = A' * transpose(sr); + mf_norm = mf_res ./ multiple_r; + echo_r_mtx(numA, :, numP) = mf_norm; + end + end + + % calculate output noise + sigma_n_o = sqrt(sigma_n^2 / multiple_r); + +end diff --git a/0 - example/wzk - RDA/rda_detection.m b/0 - example/wzk - RDA/rda_detection.m new file mode 100644 index 0000000..bb3e18c --- /dev/null +++ b/0 - example/wzk - RDA/rda_detection.m @@ -0,0 +1,47 @@ +function [ target_list_RDA, target_map ] = rda_detection(stat_RD, P_fa) +% Target detection, estimating distance, velocity, and angle intervals +% +% Usage: +% [target_list_RDA, target_map] = rda_detection(stat_RD, P_fa) +% +% Inputs: +% stat_RD: Range-Doppler statistics +% P_fa: Probability of false alarm +% +% Outputs: +% target_list_RDA: Target list include range, doppler and angle information +% target_map: Show targets in RD map + + + % parameters + [lenA, lenR, lenP] = size(stat_RD); + target_map = zeros(lenR, lenP); + + % detection threshold + d_thd = chi2inv(1 - P_fa, 2) / 2; + + % find target + target_list_RD = []; + for i = 1: lenA + RD_map = squeeze(stat_RD(i, :, :)); + detect_map = zeros(size(RD_map)); + detect_map(RD_map > d_thd) = 1; +% target_map(i, :, :) = detect_map; + [r, c] = find(detect_map); + temp_node = [r, c]; + target_list_RD = [target_list_RD; temp_node]; + end + target_list_RD_new = unique(target_list_RD, 'rows'); + + % estimate target angle interval + target_list_angle = []; + for i = 1: size(target_list_RD_new, 1) + tgt = target_list_RD_new(i, :); + stat_vct = stat_RD(:, tgt(1), tgt(2)); + [v, idx] = max(stat_vct); + target_map(tgt(1), tgt(2)) = idx; + target_list_angle = [target_list_angle; idx]; + end + target_list_RDA = [target_list_RD_new, target_list_angle]; + +end diff --git a/0 - example/wzk - RDA/spatial_FFT.m b/0 - example/wzk - RDA/spatial_FFT.m new file mode 100644 index 0000000..01ccd19 --- /dev/null +++ b/0 - example/wzk - RDA/spatial_FFT.m @@ -0,0 +1,31 @@ +function [ echo_sFFT_mtx, angle_index ] = spatial_FFT(echo_mtx, d, lambda) +% FFT for echo space domain +% +% Usage: +% [echo_sFFT_mtx, angle_index] = spatial_FFT(echo_mtx, d, lambda) +% +% Inputs: +% echo_mtx: Matrix containing the echo data +% d: Antenna spacing +% lambda: Wave length +% +% Outputs: +% echo_sFFT_mtx: Processed echo matrix +% angle_index: True value of the angle + + + % parameters + [lenA, lenR, lenP] = size(echo_mtx); + + % calculate angle + angle_index = asin(-1 * (-1/2: 1/lenA: 1/2 - 1/lenA) * lambda / d) / pi * 180; + + % spatial FFT + echo_mtx = reshape(echo_mtx, [lenA, lenR * lenP]); + echo_sFFT_mtx = zeros(size(echo_mtx)); + for i = 1: lenR * lenP + echo_sFFT_mtx(:, i) = fftshift(fft(echo_mtx(:, i))) ./ sqrt(lenA); + end + echo_sFFT_mtx = reshape(echo_sFFT_mtx, [lenA, lenR, lenP]); + +end diff --git a/0 - example/wzk - main_process/angle_est.m b/0 - example/wzk - main_process/angle_est.m new file mode 100644 index 0000000..ede069b --- /dev/null +++ b/0 - example/wzk - main_process/angle_est.m @@ -0,0 +1,43 @@ +% echo_r_mtx: range processing echo data +% target_list_RD: target list include range and doppler information +% d: antenna spacing +% fc: carrier frequency +% acc: angle estimation accuracy + +% target_list_RDA: target list include range, doppler and angle information +% RA_map: range-angle map +% THETA: scale of angle + + +function [ target_list_RDA, RA_map, THETA ] = angle_est(echo_r_mtx, target_list_RD, d, fc, acc) + + % parameters + if nargin < 5 + acc = 181; + end + [lenA, lenR, lenP] = size(echo_r_mtx); + c = 3e8; + lambda = c / fc; + + % initialization + echo_mtx = squeeze(echo_r_mtx(:, :, 1)); + + % CBF + THETA = linspace(-90, 90, round(acc)); + RA_map = zeros(length(THETA), lenR); + for i = 1: length(THETA) + a = exp((0: lenA - 1)' * -1j * 2 * pi / lambda * d * sin(THETA(i) / 180 * pi)); + RA_map(i, :) = a'* echo_mtx; + end + RA_map = transpose(RA_map); + + % calculate angle + range_angle = zeros(1, lenR); + for j = 1: lenR + [v, num] = max(abs(RA_map(j, :))); + range_angle(j) = THETA(num); + end + targets_angle = range_angle(target_list_RD(:, 2)); + target_list_RDA = [target_list_RD, transpose(targets_angle)]; + +end diff --git a/0 - example/wzk - main_process/doppler_process_CS.m b/0 - example/wzk - main_process/doppler_process_CS.m new file mode 100644 index 0000000..15f94ad --- /dev/null +++ b/0 - example/wzk - main_process/doppler_process_CS.m @@ -0,0 +1,187 @@ +% echo_r_mtx: range processing echo data +% sigma_n: input noise standard deviation +% lambda: LASSO weight +% gamma: compressed ratio +% delta: convergence normalized difference + +% echo_rd_mtx: range-doppler processing echo data +% stat_RD: range-doppler statistics + + +function [ echo_rd_mtx, stat_RD ] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta) + + % parameters + if nargin < 5 + delta = 1e-6; + end + if nargin < 4 + gamma = 0.5; + end + [lenA, lenR, M] = size(echo_r_mtx); + N = round(M / gamma); + iter_max_VAMP = 1000; + lambda_v = zeros(N, 1) + lambda; + + % generate mtx + F_ori = dftmtx(N); + F = F_ori(1:M,:); + F_inv = conj(F) / N; + + % normalization + A = (sqrt(N) * eye(M)) * F_inv; + echo_r_mtx = sqrt(N) .* echo_r_mtx; + sigma_n = sqrt(N) * sigma_n; + + % doppler matched filtering + echo_rd_mtx = zeros(lenA, lenR, N); + stat_RD = zeros(lenA, lenR, N); + for numA = 1: lenA + for numR = 1: lenR + sample = squeeze(echo_r_mtx(numA, numR, :)); + y = sample; + x_LASSO = cVAMPro(y, A, lambda_v, delta, iter_max_VAMP); + [x_d_CROD, sigma_CROD] = CROD(y, A, x_LASSO, lambda, sigma_n); + stat_RD(numA, numR, :) = abs(fftshift(x_d_CROD) / sigma_CROD).^2; + echo_rd_mtx(numA, numR, :) = fftshift(x_d_CROD); + end + end + +end + +% algorithm for LASSO +% y: measurements +% A: measurement matrix +% lambda: LASSO weight +% tau: convergence normalized difference +% Kit: maximum number of iterations +% LASSO estimator +function x_hat_wl = cVAMPro(y, A, lambda, tau, Kit) + + % Initialization + [M, N] = size(A); + gamma = M / N; + k = 0; + p = ctranspose(A) * y; + h_1 = p; + Q_1 = gamma; + tau_d = 1; + + % Iteration + while ((k < Kit) && (tau_d > tau)) + % Factorized Part + x_1 = ST(h_1, lambda, Q_1); + chi_1 = F1(x_1, lambda, Q_1); + % Message Passing + h_2 = x_1 / chi_1 - h_1; + Q_2 = 1 / chi_1 - Q_1; + % Gaussian Part + t1 = (p + h_2) / Q_2; + t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2); + x_2 = t1 - t2; + chi_2 = gamma / (Q_2 + 1) + (1 - gamma) / Q_2; + % Message Passing + h_1_next = x_2 ./ chi_2 - h_2; + Q_1_next = 1 / chi_2 - Q_2; + tau_d = norm(h_1_next - h_1, Inf) / norm(h_1_next, Inf); + k = k + 1; + % output + x_hat_wl = x_1; + % next + h_1 = h_1_next; + Q_1 = Q_1_next; + end + +end + + +% soft threshold function +% x: processing object +% thd: threshold +% y: result +function x = ST(h_1, lambda, Q_1) + [N, M] = size(h_1); + x = zeros(N, M); + + for i = 1:N + sign = h_1(i) ./ abs(h_1(i)); + diff = abs(h_1(i)) - lambda(i); + x(i) = sign .* (diff ./ Q_1) .* SF(diff); + end + +end + + +% Heaviside's step function +function v = SF(a) + + if a > 0 + v = 1; + elseif a == 0 + v = 0; % at zero points + else + v = 0; + end + +end + + +% Calculation of chi_1 +function chi_1 = F1(x_1, lambda, Q_1) + + [N, M] = size(x_1); + count = 0; + + for i = 1:N + temp = Q_1 * abs(x_1(i)) + lambda(i); + count = count + (2 - lambda(i) / temp) * SF(abs(x_1(i))); +% count = count + (2-lambda(i)/temp) * (abs(x_1(i)) > 1e-4); + end + + chi_1 = count / (2 * N * Q_1); + +end + +% calculate debiased LASSO estimator +% y: measurements +% A: measurement matrix +% x_LASSO: LASSO estimator +% lambda: LASSO weight +% sigma_n: input noise standard deviation +% x_d_CROD: debiased LASSO estimator +% sigma_CROD: equivalent noise standard deviation estimator +function [ x_d_CROD, sigma_CROD ] = CROD(y, A, x_LASSO, lambda, sigma_n) + + [m, n] = size(A); + gamma = m / n; + + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat; + + RSS = sum(abs(y - A * x_LASSO).^2)/m; + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + +end diff --git a/0 - example/wzk - main_process/doppler_process_MF.m b/0 - example/wzk - main_process/doppler_process_MF.m new file mode 100644 index 0000000..f66469b --- /dev/null +++ b/0 - example/wzk - main_process/doppler_process_MF.m @@ -0,0 +1,38 @@ +% echo_r_mtx: range processing echo data +% sigma_n: input noise standard deviation +% gamma: compressed ratio + +% echo_rd_mtx: range-doppler processing echo data +% stat_RD: range-doppler statistics + + +function [ echo_rd_mtx, stat_RD ] = doppler_process_MF(echo_r_mtx, sigma_n, gamma) + + % parameters + if nargin < 3 + gamma = 0.5; + end + [lenA, lenR, M] = size(echo_r_mtx); + N = round(M / gamma); + + % generate mtx + F_ori = dftmtx(N); + F = F_ori(1:M,:); + multiple_d = F(:,1)' * F(:,1); + + % doppler matched filtering + echo_rd_mtx = zeros(lenA, lenR, N); + for numA = 1: lenA + for numR = 1: lenR + sample = squeeze(echo_r_mtx(numA, numR, :)); + dpl_temp = transpose(F) * sample; + dpl_norm = fftshift(dpl_temp ./ multiple_d); + echo_rd_mtx(numA, numR, :) = dpl_norm; + end + end + + % calculate output noise + sigma_n_o = sqrt(sigma_n^2 / multiple_d); + stat_RD = abs(echo_rd_mtx ./ sigma_n_o).^2; + +end diff --git a/0 - example/wzk - main_process/draw_target_map.m b/0 - example/wzk - main_process/draw_target_map.m new file mode 100644 index 0000000..9fee419 --- /dev/null +++ b/0 - example/wzk - main_process/draw_target_map.m @@ -0,0 +1,21 @@ +% target_list_RDA: target list, including range, doppler and angle information +% stat_RD: range-doppler statistics + +% target_map: show targets information + + +function [ target_map ] = draw_target_map(target_list_RDA, stat_RD) + + % initialization + target_map = ones(size(stat_RD)) .* -300; + t_map = target_list_RDA(:, 1); + t_R = target_list_RDA(:, 2); + t_D = target_list_RDA(:, 3); + t_A = target_list_RDA(:, 4); + + % draw_map + for i = 1: size(target_list_RDA, 1) + target_map(t_map(i), t_R(i), t_D(i)) = t_A(i); + end + +end diff --git a/0 - example/wzk - main_process/generate_chirp_mtx.m b/0 - example/wzk - main_process/generate_chirp_mtx.m new file mode 100644 index 0000000..fc58d0f --- /dev/null +++ b/0 - example/wzk - main_process/generate_chirp_mtx.m @@ -0,0 +1,72 @@ +% PRF: pulse repetition frequency +% B: bandwidth +% fs: sampling rate +% D: duty ratio +% gamma: compression ratio + +% A: chirp matrix +% signal_t: transmitting beam + + +function [ A, signal_t ] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid) + + % parameters + if nargin < 5 + gamma = 0.5; + end + if nargin < 6 + sign_mid = 0; + end + Tr = 1 / PRF; + Tp = Tr * D; + K = B / Tp; + + % generate_signal + N = Tr * fs; + N_high = Tp * fs; + N_mtx = round(N / gamma); + + signal_t = zeros(1, N); + for i = 1: N_high + tp = i * (1 / fs) - sign_mid * N_high / fs / 2; + signal_t(1, i) = exp(1j * 2 * pi * 0.5 * K * tp .^ 2); + end + + signal_t_2fs = zeros(1, N_mtx); + for i = 1: round(N_high / gamma) + tp = (i+1) * (1 / fs / 2) - sign_mid * N_high / fs / 2; + signal_t_2fs(1, i) = exp(1j * pi * K * tp .^ 2); + end + + % generate chirp matrix + A = generate_matrix_by_signal2fs(transpose(signal_t_2fs)); + +end + + +% generate chirp matrix when gamma=0.5 +function [ mtx ] = generate_matrix_by_signal2fs(signal) + + mtx = []; + l = round(length(signal) / 2); + temp1 = signal(1:2:end); + temp2 = circshift(signal(2:2:end), 1); + if length(temp2) ~= length(temp1) + temp2 = [temp2; 0]; + end + + for i = 1:l + t1 = circshift(temp1, i-1); + t2 = circshift(temp2, i-1); + if i - 1 > 0 + t1(1:i - 1,1) = 0; + end + if i - 1 > 0 + t2(1:i - 1,1) = 0; + end + mtx = [mtx,t1,t2]; + end + +end + + diff --git a/0 - example/wzk - main_process/get_range_speed_val.m b/0 - example/wzk - main_process/get_range_speed_val.m new file mode 100644 index 0000000..3fde3fd --- /dev/null +++ b/0 - example/wzk - main_process/get_range_speed_val.m @@ -0,0 +1,29 @@ +% PRF: pulse repetition frequency +% fs: sampling rate +% fc: carrier frequency +% numP: number of pulses +% gamma_r: range compression ratio +% gamma_d: doppler compression ratio + +% distance_v: real distance value +% speed_v: real speed value + + +function [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, numP, gamma_r, gamma_d) + + % parameters + if nargin < 6 + gamma_d = 0.5; + end + if nargin < 5 + gamma_r = 0.5; + end + c = 3e8; + Tr = 1 / PRF; + lambda = c / fc; + + % calculate + distance_v = 0: (c / fs / 2 * gamma_r) : (Tr * c / 2 - c / fs / 2 * gamma_r); + speed_v = -(-PRF / 2: PRF / numP * gamma_d: PRF / 2 - PRF / numP * gamma_d) * lambda / 2 ; + +end \ No newline at end of file diff --git a/0 - example/wzk - main_process/main_cs_part.m b/0 - example/wzk - main_process/main_cs_part.m new file mode 100644 index 0000000..7a42fde --- /dev/null +++ b/0 - example/wzk - main_process/main_cs_part.m @@ -0,0 +1,84 @@ +clc +clear +close all + +%% load echo data +filename = 'Raw_Echo_60dB'; +load(['./data/', filename, '.mat']); + + +%% parameters +% ladar +PRF = 5000; +B = 5e6; +D = 0.1; +Tp = 2e-5; +fs = 5e6; +fc = 1.25e9; + +% antenna +num_antenna = 18; +d = 0.12; + +% echo +num_pulse = 64; +sigma_n = 0.1; + +% detect +P_fa = 1e-5; + +% algorithm +gamma_r = 0.5; +gamma_d = 0.5; +lambda_r = 0.005; +lambda_d = 0.15; + + +%% data processing +fprintf(['File: ', filename, '\n']); +[distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pulse, gamma_r, gamma_d); + +[A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); +fprintf(' 1: generate_chirp_mtx done.\n'); + +% first antenna +Raw_Echo_antenna1 = Raw_Echo(1, :, :); +[echo_r_mtx_cs_antenna1, sigma_n_o_cs_antenna1] = range_process_CS(Raw_Echo_antenna1, A, sigma_n, lambda_r); +fprintf(' 2: first antenna range_process_CS done.\n'); + +% first pulse +Raw_Echo_pulse1 = Raw_Echo(:, :, 1); +[echo_r_mtx_cs_pulse1, sigma_n_o_cs_pulse1] = range_process_CS(Raw_Echo_pulse1, A, sigma_n, lambda_r); +fprintf(' 3: first pulse range_process_CS done.\n'); + +[echo_rd_mtx_cs, stat_RD_cs] = doppler_process_CS(echo_r_mtx_cs_antenna1, sigma_n_o_cs_antenna1, lambda_d, gamma_d); +fprintf(' 4: doppler_process_CS done.\n'); + +[target_list_RD_cs, target_map_cs] = rd_detection(stat_RD_cs, P_fa); +fprintf(' 5: rd_detection done.\n'); + +[target_list_RDA, RA_map] = angle_est(echo_r_mtx_cs_pulse1, target_list_RD_cs, d, fc); +fprintf(' 6: angle_est done.\n'); + + +%% show result +element = 1; +list_r_idx = find(target_list_RDA(:, 1) == element); + +target_list_RDA(:, 2) = distance_v(target_list_RDA(:, 2)); +target_list_RDA(:, 3) = speed_v(target_list_RDA(:, 3)); +figure(1) +scatter3(target_list_RDA(list_r_idx,2), target_list_RDA(list_r_idx,3),... + target_list_RDA(list_r_idx,4), 'filled', 'o') + +% xlim([17000, 19000]) +% ylim([80, 160]) +% zlim([-90, 90]) + + + +%% save +echo_rd_mtx = echo_rd_mtx_cs; +save(['./output/', filename, '_cs_part.mat'],... + 'target_list_RDA', 'echo_rd_mtx', 'RA_map', 'P_fa'); + diff --git a/0 - example/wzk - main_process/main_mf.m b/0 - example/wzk - main_process/main_mf.m new file mode 100644 index 0000000..dd19dcf --- /dev/null +++ b/0 - example/wzk - main_process/main_mf.m @@ -0,0 +1,67 @@ +clc +clear +close all + +%% load echo data +filename = 'Raw_Echo_60dB'; +load(['./data/', filename, '.mat']); + + +%% parameters +% ladar +PRF = 5000; +B = 5e6; +D = 0.1; +Tp = 2e-5; +fs = 5e6; +fc = 1.25e9; + +% antenna +num_antenna = 18; +d = 0.12; + +% echo +num_pulse = 64; +sigma_n = 0.1; + +% detect +P_fa = 1e-4; + +% algorithm +gamma_r = 0.5; +gamma_d = 0.5; + + +%% data processing +[distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pulse, gamma_r, gamma_d); + +[A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + +[echo_r_mtx_mf, sigma_n_o_mf] = range_process_MF(Raw_Echo, A, sigma_n); + +[echo_rd_mtx_mf, stat_RD_mf] = doppler_process_MF(echo_r_mtx_mf, sigma_n_o_mf, gamma_d); + +[target_list_RD_mf, target_map_mf] = rd_detection(stat_RD_mf, P_fa); + +[target_list_RDA, RA_map] = angle_est(echo_r_mtx_mf, target_list_RD_mf, d, fc); + + +%% show result +element = 1; +list_r_idx = find(target_list_RDA(:, 1) == element); + +target_list_RDA(:, 2) = distance_v(target_list_RDA(:, 2)); +target_list_RDA(:, 3) = speed_v(target_list_RDA(:, 3)); +figure(1) +scatter3(target_list_RDA(list_r_idx,2), target_list_RDA(list_r_idx,3),... + target_list_RDA(list_r_idx,4), 'filled', 'o') +% xlim([17000, 19000]) +% ylim([80, 160]) +% zlim([-90, 90]) + + +%% save +echo_rd_mtx = echo_rd_mtx_mf; +save(['./output/', filename, '_mf.mat'],... + 'target_list_RDA', 'echo_rd_mtx', 'RA_map', 'P_fa'); + diff --git a/0 - example/wzk - main_process/main_mf_part.m b/0 - example/wzk - main_process/main_mf_part.m new file mode 100644 index 0000000..bb994ac --- /dev/null +++ b/0 - example/wzk - main_process/main_mf_part.m @@ -0,0 +1,73 @@ +clc +clear +close all + +%% load echo data +filename = 'Raw_Echo_60dB'; +load(['./data/', filename, '.mat']); + + +%% parameters +% ladar +PRF = 5000; +B = 5e6; +D = 0.1; +Tp = 2e-5; +fs = 5e6; +fc = 1.25e9; + +% antenna +num_antenna = 18; +d = 0.12; + +% echo +num_pulse = 64; +sigma_n = 0.1; + +% detect +P_fa = 1e-4; + +% algorithm +gamma_r = 0.5; +gamma_d = 0.5; + + +%% data processing +[distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pulse, gamma_r, gamma_d); + +[A, signal_t] = generate_chirp_mtx(PRF, B, fs, D); + +% first antenna +Raw_Echo_antenna1 = Raw_Echo(1, :, :); +[echo_r_mtx_mf_antenna1, sigma_n_o_mf_antenna1] = range_process_MF(Raw_Echo_antenna1, A, sigma_n); + +% first pulse +Raw_Echo_pulse1 = Raw_Echo(:, :, 1); +[echo_r_mtx_mf_pulse1, sigma_n_o_mf_pulse1] = range_process_MF(Raw_Echo_pulse1, A, sigma_n); + +[echo_rd_mtx_mf, stat_RD_mf] = doppler_process_MF(echo_r_mtx_mf_antenna1, sigma_n_o_mf_antenna1, gamma_d); + +[target_list_RD_mf, target_map_mf] = rd_detection(stat_RD_mf, P_fa); + +[target_list_RDA, RA_map] = angle_est(echo_r_mtx_mf_pulse1, target_list_RD_mf, d, fc); + + +%% show result +element = 1; +list_r_idx = find(target_list_RDA(:, 1) == element); + +target_list_RDA(:, 2) = distance_v(target_list_RDA(:, 2)); +target_list_RDA(:, 3) = speed_v(target_list_RDA(:, 3)); +figure(2) +scatter3(target_list_RDA(list_r_idx,2), target_list_RDA(list_r_idx,3),... + target_list_RDA(list_r_idx,4), 'filled', 'o') +% xlim([17000, 19000]) +% ylim([80, 160]) +% zlim([-90, 90]) + + +%% save +echo_rd_mtx = echo_rd_mtx_mf; +save(['./output/', filename, '_mf_part.mat'],... + 'target_list_RDA', 'echo_rd_mtx', 'RA_map', 'P_fa'); + diff --git a/0 - example/wzk - main_process/plot_result.m b/0 - example/wzk - main_process/plot_result.m new file mode 100644 index 0000000..9d5b7aa --- /dev/null +++ b/0 - example/wzk - main_process/plot_result.m @@ -0,0 +1,47 @@ +clc +clear +close all + + +Fontsize = 18; +plot_width = 800; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + + +Pfa_set = '1e-4'; +method = 'cs'; +if_part = 1; +SNR = 10; + +if if_part == 1 + filename = ['Raw_Echo_', num2str(SNR), 'dB_', method, '_part']; +else + filename = ['Raw_Echo_', num2str(SNR), 'dB_', method]; +end + +load(['./output/Pfa', Pfa_set, '/', filename, '.mat']); + +%% +element = 1; +list_r_idx = find(target_list_RDA(:, 1) == element); + +figure(1) +scatter3(target_list_RDA(list_r_idx,2), target_list_RDA(list_r_idx,3),... + target_list_RDA(list_r_idx,4), 'filled', 'o') +xlabel('Range(m)'); +ylabel('Speed(m/s)'); +zlabel('Angle(°)'); +zlim([-90, 90]); +ylim([80, 160]); +xlim([17000, 21000]); +set(gca, 'FontSize', Fontsize); +title(['Target detected by ', upper(method), ' under Pfa=', Pfa_set]) +set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]); +set(gca,'fontsize',18,'fontname','Times'); + + + + + diff --git a/0 - example/wzk - main_process/range_process_CS.m b/0 - example/wzk - main_process/range_process_CS.m new file mode 100644 index 0000000..865cdcd --- /dev/null +++ b/0 - example/wzk - main_process/range_process_CS.m @@ -0,0 +1,146 @@ +% echo_mtx: echo data +% A: chirp matrix +% sigma_n: input noise standard deviation +% lambda: LASSO weight +% delta: convergence normalized difference + +% echo_r_mtx: range processing echo data +% signal_n_o: output noise standard deviation + + +function [ echo_r_mtx, sigma_n_o ] = range_process_CS(echo_mtx, A, sigma_n, lambda, delta) + + % parameters + if nargin < 5 + delta = 2e-7; + end + [M, N] = size(A); + [lenA, M, lenP] = size(echo_mtx); + + % normalization + J1 = A*A'; + lambda_J=eig(J1); + A = A / sqrt(lambda_J(end)); + echo_mtx = echo_mtx ./ sqrt(lambda_J(end)); + sigma_n = sigma_n / sqrt(lambda_J(end)); + + % compressed sensing + echo_r_mtx = zeros(lenA, N, lenP); + for numA = 1: lenA + for numP = 1: lenP + sr = echo_mtx(numA, :, numP); + y = transpose(sr); + % LASSO + x_FISTA = FISTA(y, A, lambda, delta); + % debiased LASSO + [x_d, sigma_w] = cal_debiased_LASSO(x_FISTA, A, y, lambda, sigma_n); + echo_r_mtx(numA, :, numP) = x_d; + end + end + + % calculate output noise + sigma_n_o = sigma_w; + +end + + +% algorithm for LASSO +% y: measurements +% A: measurement matrix +% lambda: LASSO weight +% delta: convergence normalized difference +% z: LASSO estimator +function [z] = FISTA(y, A, lambda, delta) + + x_pre = A'*y; + t = 1; + z = x_pre; + z_pre = z; + t_pre = t; + N = size(A, 2); + diff = 1; + E = eig(A'*A); + L = E(end); + temp1 = A'*y/L; + temp2 = eye(N) - A'*A/L; + k = 0; + + while((diff > delta) && (k < 1000)) + temp = temp1 + temp2 * z_pre; + x = sft_thd(temp, lambda/L); + t = 0.5*(1 + sqrt(1+4*t_pre*t_pre)); + z = x + (x - x_pre) * (t_pre-1) / t; + diff = mean(abs(z_pre - z)); + x_pre = x; + z_pre = z; + t_pre = t; + k = k + 1; + end + +end + + +% soft threshold function +% x: processing object +% thd: threshold +% y: result +function y = sft_thd(x, thd) + + if isequal(size(x), size(thd)) + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd); + else + tmp = abs(x); + y = x; + y(tmp <= thd) = 0; + y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd); + end + +end + + +% calculate debiased LASSO estimator +% x: LASSO estimator +% A: measurement matrix +% y: measurements +% lambda: LASSO weight +% sigma_n: input noise standard deviation +% x_d: debiased LASSO estimator +% sigma_d: equivalent noise standard deviation +function [x_d, sigma_d] = cal_debiased_LASSO(x, A, y, lambda, sigma) + + [M, N] = size(A); + gamma = M/N; + hat_Q1 = gamma; + [~, D] = eig(A'*A); + d = diag(D); + + diff = 1; + T = 1000; + t = 0; + + while (t < T) && (diff > 1e-6) + Q1_pre = hat_Q1; + rho = mean((2 - lambda./(hat_Q1*abs(x) + lambda)).*(abs(x) > 1e-4))/2; + hat_Q1 = rho/mean(1./(d + (1-rho)*hat_Q1/rho)); + diff = abs(Q1_pre - hat_Q1); + t = t+1; + end + + x_d = x + 1/hat_Q1*A'*(y - A*x); + + chi = rho/hat_Q1; + hat_Q2 = 1/chi - hat_Q1; + + t = -hat_Q2; + t_prime = -1/mean((1./(d+hat_Q2)).^2); + G_prime = t + 1/chi; + G_wprime = t_prime + 1/chi/chi; + RSS = sum(abs(y - A*x).^2)/M; + hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +... + (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2; + sigma_d = sqrt(2*hat_chi)/hat_Q1; + +end diff --git a/0 - example/wzk - main_process/range_process_MF.m b/0 - example/wzk - main_process/range_process_MF.m new file mode 100644 index 0000000..977ceed --- /dev/null +++ b/0 - example/wzk - main_process/range_process_MF.m @@ -0,0 +1,30 @@ +% echo_mtx: echo data +% A: chirp matrix +% sigma_n: input noise standard deviation + +% echo_r_mtx: range processing echo data +% signal_n_o: output noise standard deviation + + +function [ echo_r_mtx, sigma_n_o ] = range_process_MF(echo_mtx, A, sigma_n) + + % parameters + multiple_r = A(:, 1)' * A(:, 1); + [M, N] = size(A); + [lenA, M, lenP] = size(echo_mtx); + + % matched filtering + echo_r_mtx = zeros(lenA, N, lenP); + for numA = 1: lenA + for numP = 1: lenP + sr = echo_mtx(numA, :, numP); + mf_res = A' * transpose(sr); + mf_norm = mf_res ./ multiple_r; + echo_r_mtx(numA, :, numP) = mf_norm; + end + end + + % calculate output noise + sigma_n_o = sqrt(sigma_n^2 / multiple_r); + +end diff --git a/0 - example/wzk - main_process/rd_detection.m b/0 - example/wzk - main_process/rd_detection.m new file mode 100644 index 0000000..eb9c4f3 --- /dev/null +++ b/0 - example/wzk - main_process/rd_detection.m @@ -0,0 +1,32 @@ +% stat_RD: range-doppler statistics +% P_fa: false alarm rate + +% target_map: show targets in RD map +% target list include range and doppler information + + +function [ target_list_RD, target_map ] = rd_detection(stat_RD, P_fa) + + % parameters + target_map = zeros(size(stat_RD)); + [lenA, lenR, lenP] = size(stat_RD); + + % detection threshold + d_thd = chi2inv(1 - P_fa, 2) / 2; + + % find target + target_list_RD = []; + for i = 1: lenA + RD_map = squeeze(stat_RD(i, :, :)); + detect_map = zeros(size(RD_map)); + detect_map(RD_map > d_thd) = 1; + target_map(i, :, :) = detect_map; + [r, c] = find(detect_map); + temp_list = zeros(length(r), 3); + temp_list(:, 1) = i; + temp_list(:, 2) = r; + temp_list(:, 3) = c; + target_list_RD = [target_list_RD; temp_list]; + end + +end diff --git a/0 - example/wzk - 目标检测/generate_chirp_mtx.m b/0 - example/wzk - 目标检测/generate_chirp_mtx.m new file mode 100644 index 0000000..757cb92 --- /dev/null +++ b/0 - example/wzk - 目标检测/generate_chirp_mtx.m @@ -0,0 +1,71 @@ +% PRF: pulse repetition frequency +% B: bandwidth +% fs: sampling rate +% D: duty ratio +% gamma: compression ratio + +% A: chirp matrix +% signal_t: transmitting beam + + +function [ A, signal_t ] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid) + + % parameters + if nargin < 5 + gamma = 0.5; + end + if nargin < 6 + sign_mid = 0; + end + Tr = 1 / PRF; + Tp = Tr * D; + K = B / Tp; + + % generate_signal + N = Tr * fs; + N_high = Tp * fs; + N_mtx = round(N / gamma); + + signal_t = zeros(1, N); + for i = 1: N_high + tp = i * (1 / fs) - sign_mid * N_high / fs / 2; + signal_t(1, i) = exp(1j * 2 * pi * 0.5 * K * tp .^ 2); + end + + signal_t_2fs = zeros(1, N_mtx); + for i = 1: round(N_high / gamma) + tp = (i+1) * (1 / fs / 2) - sign_mid * N_high / fs / 2; + signal_t_2fs(1, i) = exp(1j * pi * K * tp .^ 2); + end + + % generate chirp matrix + A = generate_matrix_by_signal2fs(transpose(signal_t_2fs)); + +end + + +function [ mtx ] = generate_matrix_by_signal2fs(signal) + + mtx = []; + l = round(length(signal) / 2); + temp1 = signal(1:2:end); + temp2 = circshift(signal(2:2:end), 1); + if length(temp2) ~= length(temp1) + temp2 = [temp2; 0]; + end + + for i = 1:l + t1 = circshift(temp1, i-1); + t2 = circshift(temp2, i-1); + if i - 1 > 0 + t1(1:i - 1,1) = 0; + end + if i - 1 > 0 + t2(1:i - 1,1) = 0; + end + mtx = [mtx,t1,t2]; + end + +end + + diff --git a/0 - example/wzk - 目标检测/node_detect.m b/0 - example/wzk - 目标检测/node_detect.m new file mode 100644 index 0000000..ef6c326 --- /dev/null +++ b/0 - example/wzk - 目标检测/node_detect.m @@ -0,0 +1,87 @@ +clc +clear +close all; + +SNR = 10; +sigma_n = 0.1; +P_fa = [1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1]; +len_P_fa = length(P_fa); + +h0_rep_time = 99 * 1e4; +h1_rep_time = 1e5; +% h0_rep_time = 1000; +% h1_rep_time = 1000; + +a = sqrt(10^(SNR/10) * sigma_n^2); + +P_fa_node_cnt = zeros(len_P_fa, h0_rep_time); + +parfor rep = 1: h0_rep_time +% for rep = 1: h0_rep_time + + %% 噪声处理 + noise = random('Normal', 0, sigma_n/sqrt(2), 1, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, 1); + + y = noise; + + stat = abs(y)^2; + + kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2; + + for cnt_h_th = 1: len_P_fa + P_fa_node_cnt(cnt_h_th, rep) = stat > kd(cnt_h_th); + end + + fprintf('Pfa-%d\n', rep); + +end + +P_fa_node = mean(P_fa_node_cnt, 2); + +figure(1) +loglog(P_fa,P_fa_node) +title('Pfa') + + + +P_d_node_cnt = zeros(len_P_fa, h1_rep_time); + +parfor rep = 1: h1_rep_time +% for rep = 1: h1_rep_time + + %% 噪声处理 + noise = random('Normal', 0, sigma_n/sqrt(2), 1, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, 1); + + y = a + noise; + + stat = abs(y)^2; + + kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2; + + for cnt_h_th = 1: len_P_fa + P_d_node_cnt(cnt_h_th, rep) = stat > kd(cnt_h_th); + end + + fprintf('Pd-%d\n', rep); + +end + +P_d_node = mean(P_d_node_cnt, 2); + +figure(2) +loglog(P_fa,P_d_node) +title('Pd') + +figure(3) +loglog(P_fa_node,P_d_node) +title('ROC') + +save node_detect.mat ... + SNR... + P_fa... + P_fa_node... + P_d_node... + sigma_n... + h0_rep_time... + h1_rep_time... + a; \ No newline at end of file diff --git a/0 - example/wzk - 目标检测/node_detect_norm.m b/0 - example/wzk - 目标检测/node_detect_norm.m new file mode 100644 index 0000000..ef6c326 --- /dev/null +++ b/0 - example/wzk - 目标检测/node_detect_norm.m @@ -0,0 +1,87 @@ +clc +clear +close all; + +SNR = 10; +sigma_n = 0.1; +P_fa = [1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1]; +len_P_fa = length(P_fa); + +h0_rep_time = 99 * 1e4; +h1_rep_time = 1e5; +% h0_rep_time = 1000; +% h1_rep_time = 1000; + +a = sqrt(10^(SNR/10) * sigma_n^2); + +P_fa_node_cnt = zeros(len_P_fa, h0_rep_time); + +parfor rep = 1: h0_rep_time +% for rep = 1: h0_rep_time + + %% 噪声处理 + noise = random('Normal', 0, sigma_n/sqrt(2), 1, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, 1); + + y = noise; + + stat = abs(y)^2; + + kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2; + + for cnt_h_th = 1: len_P_fa + P_fa_node_cnt(cnt_h_th, rep) = stat > kd(cnt_h_th); + end + + fprintf('Pfa-%d\n', rep); + +end + +P_fa_node = mean(P_fa_node_cnt, 2); + +figure(1) +loglog(P_fa,P_fa_node) +title('Pfa') + + + +P_d_node_cnt = zeros(len_P_fa, h1_rep_time); + +parfor rep = 1: h1_rep_time +% for rep = 1: h1_rep_time + + %% 噪声处理 + noise = random('Normal', 0, sigma_n/sqrt(2), 1, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, 1); + + y = a + noise; + + stat = abs(y)^2; + + kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2; + + for cnt_h_th = 1: len_P_fa + P_d_node_cnt(cnt_h_th, rep) = stat > kd(cnt_h_th); + end + + fprintf('Pd-%d\n', rep); + +end + +P_d_node = mean(P_d_node_cnt, 2); + +figure(2) +loglog(P_fa,P_d_node) +title('Pd') + +figure(3) +loglog(P_fa_node,P_d_node) +title('ROC') + +save node_detect.mat ... + SNR... + P_fa... + P_fa_node... + P_d_node... + sigma_n... + h0_rep_time... + h1_rep_time... + a; \ No newline at end of file diff --git a/0 - example/wzk - 目标检测/plot_node_detect.m b/0 - example/wzk - 目标检测/plot_node_detect.m new file mode 100644 index 0000000..7623f44 --- /dev/null +++ b/0 - example/wzk - 目标检测/plot_node_detect.m @@ -0,0 +1,42 @@ +clear; +close all; +clc; + +load node_detect.mat; + +Fontsize = 18; +plot_width = 850; +plot_height = 600; +Linewidth = 2; +Markersize = 8; + + + +%% plot +figure(1); +loglog(P_fa, P_fa_node, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +grid on; +legend('SNR = 10dB'); +xlabel('P_{fa} Set'); +ylabel('Actual P_{fa}'); +ylim([1e-5,1]); +set(gca, 'FontSize', Fontsize); +title('P_{fa}') +set(gcf, 'position', [200, 300, plot_width, plot_height]); +set(gca,'fontsize',20,'fontname','Times'); + +figure(3); +semilogx(P_fa_node, P_d_node, '-o', ... + 'Linewidth', Linewidth, ... + 'MarkerSize', Markersize); +grid on; +legend('SNR = 10dB'); +xlabel('Actual P_{fa}'); +ylabel('P_{d}'); +xlim([9.99e-6,1]); +set(gca, 'FontSize', Fontsize); +title('ROC') +set(gcf, 'position', [200, 300, plot_width, plot_height]); +set(gca,'fontsize',20,'fontname','Times'); \ No newline at end of file