檢測(cè)原理與MATLAB仿真:從Q函數(shù)到蒙特卡洛性能分析)
簡(jiǎn)介這份MATLAB源碼包圍繞參數(shù)已知條件下的廣義似然比檢驗(yàn)GLRT信號(hào)檢測(cè)問題編寫適合信號(hào)處理、統(tǒng)計(jì)檢測(cè)課程學(xué)習(xí)者與算法入門者參考。資源共4個(gè)文件壓縮包僅7KB包含3個(gè)m腳本和1張結(jié)果圖drawpd.m用于繪制兩種假設(shè)下的概率密度函數(shù)Q.m與Qinv.m涉及Q函數(shù)及逆Q函數(shù)可用于計(jì)算檢測(cè)閾值和誤警概率結(jié)果.bmp直觀展示了仿真輸出的檢測(cè)性能。目前已有1107人學(xué)習(xí)下載。通過研讀和運(yùn)行源碼可系統(tǒng)掌握GLRT從模型建立、對(duì)數(shù)似然比構(gòu)造到閾值判決的完整流程同時(shí)學(xué)習(xí)如何在MATLAB中實(shí)現(xiàn)統(tǒng)計(jì)檢測(cè)仿真、繪制PDF曲線并評(píng)估檢測(cè)效果。該小體積實(shí)例結(jié)構(gòu)清晰、便于調(diào)試適合作為課程實(shí)驗(yàn)或自學(xué)GLRT算法的起步模板。1. 為什么參數(shù)已知時(shí) GLRT 是信號(hào)檢測(cè)的首選工具做信號(hào)檢測(cè)的人大概都遇到過這樣的場(chǎng)景噪聲統(tǒng)計(jì)特性已知信號(hào)幅度、相位、到達(dá)時(shí)間也由先驗(yàn)信息給定唯一不確定的就是“此刻到底有沒有信號(hào)”。這時(shí)如果還去套能量檢測(cè)或者匹配濾波的固定判決門限往往會(huì)在低信噪比下?lián)p失檢測(cè)性能而 GLRT廣義似然比檢驗(yàn)恰好能在這個(gè)前提下把檢測(cè)問題轉(zhuǎn)化為一個(gè)可解析的似然比比較過程。它不要求對(duì)未知參數(shù)做積分或貝葉斯平均而是用最大似然估計(jì)值代入似然比從而把復(fù)合假設(shè)檢驗(yàn)簡(jiǎn)化成單邊閾值判決。對(duì)做雷達(dá)目標(biāo)檢測(cè)、通信同步頭捕獲、生物醫(yī)學(xué)信號(hào)事件判別的工程師來說這套 MATLAB 仿真代碼的價(jià)值在于它把從假設(shè)建模、Q 函數(shù)查閾值到蒙特卡洛畫曲線的完整鏈路都放在你面前既能驗(yàn)證理論曲線又能直接改成自己的數(shù)據(jù)格式。2. GLRT 檢測(cè)原理與假設(shè)檢驗(yàn)建模2.1 從 Neyman-Pearson 到 GLRT 的演進(jìn)經(jīng)典檢測(cè)理論里當(dāng)噪聲是零均值高斯白噪聲、信號(hào)波形完全已知時(shí)最優(yōu)檢測(cè)器是匹配濾波加固定閾值也就是 Neyman-Pearson 準(zhǔn)則下的似然比檢驗(yàn)。似然比定義為兩個(gè)假設(shè)下觀測(cè)數(shù)據(jù)概率密度函數(shù)的比值L(x) p(x|H1) / p(x|H0)當(dāng) L(x) 超過某個(gè)閾值時(shí)判 H1否則判 H0。此時(shí)閾值由虛警概率決定而虛警概率又依賴噪聲方差和信號(hào)能量計(jì)算并不復(fù)雜。但實(shí)際工程里“完全已知”很少成立。更多情況是信號(hào)里還藏著幾個(gè)未知參數(shù)比如直流偏置、正弦波的初相、脈沖信號(hào)的到達(dá)時(shí)刻。如果把這些參數(shù)當(dāng)作隨機(jī)變量去積分就是貝葉斯檢測(cè)如果不知道它們的先驗(yàn)分布GLRT 就成了務(wù)實(shí)的選擇先用最大似然估計(jì)把未知參數(shù)估出來再代入似然比。從原理上看GLRT 是給每個(gè)可能的參數(shù)值做一次匹配再取最大輸出和門限比較性能雖然略遜于參數(shù)已知的最優(yōu)檢測(cè)但在未知參數(shù)維度不高時(shí)損失極小。參數(shù)全部已知的 GLRT 是這個(gè)框架的退化情形。這時(shí)最大似然估計(jì)退化成已知值GLRT 統(tǒng)計(jì)量就等價(jià)于標(biāo)準(zhǔn)似然比。但這個(gè)仿真項(xiàng)目特意保留了 Q.m、Qinv.m 這兩個(gè)模塊說明它不是簡(jiǎn)單調(diào)一次normcdf就完事而是把閾值計(jì)算、逆高斯分位數(shù)求解也拆成了獨(dú)立函數(shù)方便以后往參數(shù)未知的方向擴(kuò)展。2.2 參數(shù)已知情形下的似然比構(gòu)造假設(shè)觀測(cè)向量是 N 維的H0 下 x wH1 下 x s w其中 w 是零均值協(xié)方差矩陣為 σ2I 的高斯噪聲s 是已知信號(hào)向量。兩個(gè)假設(shè)下的概率密度分別為p(x|H0) (1/(2πσ2)^(N/2)) * exp(-||x||2/(2σ2)) p(x|H1) (1/(2πσ2)^(N/2)) * exp(-||x-s||2/(2σ2))兩邊取對(duì)數(shù)再相減得到對(duì)數(shù)似然比ln L(x) (1/(2σ2)) * (||x||2 - ||x-s||2) (1/σ2) * (x?s - ||s||2/2)去掉與數(shù)據(jù)無(wú)關(guān)的常數(shù)項(xiàng)之后檢測(cè)統(tǒng)計(jì)量就變成T(x) x?s也就是觀測(cè)與已知信號(hào)的內(nèi)積。這個(gè)結(jié)果和匹配濾波器是一致的。GLRT 在這里的“廣義”體現(xiàn)在如果 s 里某些參數(shù)未知那么 T(x) 里會(huì)用這些參數(shù)的最大似然估計(jì)去替換從而得到T_glrt(x) max_θ x?s(θ)。在 MATLAB 里構(gòu)造這個(gè)統(tǒng)計(jì)量不需要循環(huán)。假設(shè)噪聲方差已知為sigma2信號(hào)向量為s觀測(cè)量為x一行代碼就能算出統(tǒng)計(jì)量T (x(:) * s(:)) / sqrt(sigma2 * (s(:) * s(:))); % 歸一化相關(guān)系數(shù)形式歸一化之后T 在 H0 下服從標(biāo)準(zhǔn)正態(tài)分布在 H1 下服從均值為sqrt(SNR)的正態(tài)分布。這樣閾值可以直接從標(biāo)準(zhǔn)正態(tài)分位數(shù)得到也就是 Q 函數(shù)的逆。這也是項(xiàng)目里出現(xiàn) Qinv.m 的核心原因。參數(shù)說明x(:)與s(:)都強(qiáng)制拉成列向量避免行向量轉(zhuǎn)置錯(cuò)誤除以信號(hào)能量開根號(hào)相當(dāng)于把匹配濾波器歸一化讓統(tǒng)計(jì)量的方差恒為 1門限只和虛警率有關(guān)不再依賴信號(hào)幅度。2.3 閾值確定與 Q 函數(shù)的作用檢測(cè)門限由虛警概率Pfa決定。定義 Q 函數(shù)為標(biāo)準(zhǔn)正態(tài)分布右尾概率Q(z) ∫_z^∞ (1/√(2π)) exp(-t2/2) dt如果希望虛警率不超過Pfa則門限應(yīng)當(dāng)滿足λ Q?1(Pfa)。在 MATLAB 里erfc函數(shù)和 Q 函數(shù)的關(guān)系是Q(z) 0.5 * erfc(z/√2)所以逆函數(shù)可以寫成lambda sqrt(2) * erfcinv(2 * Pfa); % 由虛警率反推判決門限項(xiàng)目里的Qinv.m很可能就是封裝了這一行變換。這樣做的優(yōu)勢(shì)是門限獨(dú)立于信號(hào)波形只要信號(hào)能量歸一化做好同一套門限可以復(fù)用到任意波形檢測(cè)。要注意的是如果噪聲方差未知或者不是白噪聲那么門限計(jì)算里還要引入噪聲協(xié)方差矩陣的 Cholesky 分解做白化這是后續(xù)擴(kuò)展的方向之一。下表總結(jié)了不同虛警概率下對(duì)應(yīng)的門限值和等效 SNR 檢測(cè)門限方便快速核對(duì)仿真參數(shù)是否設(shè)置合理PfaQinv(Pfa) 門限 λ所需 SNR(dB) 約 90% 檢測(cè)率0.012.3263約 4.3 dB0.0013.0902約 6.0 dB1e-54.2649約 8.5 dB1e-75.1993約 10.4 dB這里的“所需 SNR”是在單次采樣檢測(cè)場(chǎng)景下把檢測(cè)概率公式Pd Q(λ - sqrt(SNR))反推得到的。仿真時(shí)如果發(fā)現(xiàn)檢測(cè)曲線和理論值偏差超過 0.5 dB首先要檢查門限是否用了未歸一化的統(tǒng)計(jì)量其次檢查蒙特卡洛次數(shù)是否足夠。3. MATLAB 仿真實(shí)現(xiàn)drawpd.m、Q.m 與 Qinv.m 拆解3.1 仿真框架與文件分工這套源碼的文件結(jié)構(gòu)非常典型適合作為檢測(cè)仿真的骨架。Q.m和Qinv.m是數(shù)值函數(shù)庫(kù)專門處理高斯 Q 函數(shù)及其逆運(yùn)算drawpd.m負(fù)責(zé)繪制概率密度曲線用來直觀展示 H0 和 H1 下統(tǒng)計(jì)量的分布結(jié)果.bmp是運(yùn)行后的輸出圖一般包含檢測(cè)概率曲線或 PDF 對(duì)比圖。這種分工方式把一個(gè)完整的檢測(cè)仿真拆成了“統(tǒng)計(jì)工具”和“場(chǎng)景腳本”后續(xù)做其他檢測(cè)器時(shí)可以直接復(fù)用Q.m和Qinv.m。從工程角度看把 Q 函數(shù)單獨(dú)抽出來是個(gè)好習(xí)慣。MATLAB 自帶的normcdf或者erfc雖然也能用但erfc對(duì)大數(shù)值參數(shù)容易下溢而自實(shí)現(xiàn)時(shí)通常會(huì)做換元或?qū)?shù)域處理。實(shí)際項(xiàng)目里我一般會(huì)再包一層讓Qinv支持向量輸入方便一次批量計(jì)算多個(gè)門限。3.2 Q 函數(shù)及逆函數(shù)的數(shù)值實(shí)現(xiàn)標(biāo)準(zhǔn)正態(tài)分布的 Q 函數(shù)可以用補(bǔ)誤差函數(shù)表達(dá)但為了數(shù)值穩(wěn)定性我建議對(duì)極限情況做保護(hù)。一個(gè)實(shí)用的Q.m實(shí)現(xiàn)如下function y Q(x) % 標(biāo)準(zhǔn)正態(tài)分布右尾概率即 P(Z x) % 輸入 x 可以是標(biāo)量或向量 % 內(nèi)部利用 erfc 精度高且支持向量化 x x(:); y 0.5 * erfc(x / sqrt(2)); % 對(duì)極端值做飽和處理避免向上/向下溢出產(chǎn)生 NaN y(x 8) 0; y(x -8) 1; endQinv.m則是對(duì)上式的逆運(yùn)算輸入概率值輸出對(duì)應(yīng)分位數(shù)function x Qinv(p) % Q 函數(shù)的逆給定右尾概率 p返回分位數(shù) x滿足 Q(x) p % 輸入 p 應(yīng)該在 (0,1) 區(qū)間內(nèi) p(p 0 | p 1) NaN; % 非法概率直接置 NaN便于調(diào)用方發(fā)現(xiàn)錯(cuò)誤 x sqrt(2) * erfcinv(2 * p); end邏輯說明由于erfc在參數(shù)絕對(duì)值較大時(shí)可能出現(xiàn)數(shù)值飽和而實(shí)際仿真中門限通常不會(huì)超過 6兩個(gè)函數(shù)的保護(hù)分支大多數(shù)情況下不會(huì)觸發(fā)。但留著它可以避免蒙特卡洛循環(huán)里因?yàn)閭€(gè)別異常樣本導(dǎo)致整個(gè)仿真崩潰。參數(shù)說明p建議取值范圍在 1e-7 到 0.5 之間如果虛警率低于 1e-7建議改用對(duì)數(shù)域遞推因?yàn)榇藭r(shí)erfcinv的參數(shù)接近 2浮點(diǎn)分辨率可能不夠。3.3 drawpd.m 中的概率密度繪制與可視化drawpd.m的核心作用是畫出 H0 和 H1 兩個(gè)假設(shè)下檢測(cè)統(tǒng)計(jì)量的理論 PDF讓使用者先看清楚“兩個(gè)分布的重疊程度”再跑蒙特卡洛。假設(shè)我們用的是歸一化內(nèi)積統(tǒng)計(jì)量那么 H0 下 T ~ N(0,1)H1 下 T ~ N(√SNR, 1)。理論 PDF 可以直接用normpdf計(jì)算function drawpd(SNR_dB) % 繪制 H0 和 H1 下的統(tǒng)計(jì)量概率密度曲線 % SNR_dB信噪比單位 dB snr 10^(SNR_dB/10); t -4:0.01:8; % 統(tǒng)計(jì)量取值范圍 pdf0 normpdf(t, 0, 1); % H0均值0方差1 pdf1 normpdf(t, sqrt(snr), 1); % H1均值 sqrt(SNR)方差1 plot(t, pdf0, b-, LineWidth, 1.5); hold on; plot(t, pdf1, r-, LineWidth, 1.5); xlabel(檢測(cè)統(tǒng)計(jì)量 T); ylabel(概率密度); legend(H0: 純?cè)肼? H1: 信號(hào)噪聲); title([PDF對(duì)比, SNR , num2str(SNR_dB), dB]); grid on; end這段代碼的關(guān)鍵參數(shù)是統(tǒng)計(jì)量的均值。sqrt(snr)來自歸一化內(nèi)積的期望推導(dǎo)當(dāng)信號(hào)能量為Es、噪聲方差為σ2歸一化后的 SNR 就是Es/σ2開根號(hào)后即為 H1 下統(tǒng)計(jì)量的均值偏移量。用-4:0.01:8作為橫軸范圍是因?yàn)殚T限通常在 2~5 之間且 H1 的均值在低 SNR 時(shí)仍在 2 附近橫軸范圍留夠余量才能看出右尾差異。drawpd.m里的繪圖語(yǔ)句不必和這個(gè)完全一致但核心邏輯一定是“把兩個(gè)假設(shè)的分布畫在同一坐標(biāo)系同時(shí)標(biāo)出判決門限”否則 PDF 圖對(duì)參數(shù)選擇的指導(dǎo)意義會(huì)大打折扣。3.4 主流程代碼與參數(shù)設(shè)置把上面的模塊串起來一個(gè)完整的參數(shù)已知 GLRT 檢測(cè)仿真主流程如下%% 參數(shù)初始化 N 10; % 采樣點(diǎn)數(shù) Pfa 1e-3; % 期望虛警率 MC 100000; % 蒙特卡洛次數(shù) SNR_dB 0:2:10; % 仿真信噪比范圍 % 已知信號(hào)單位能量矩形脈沖 s ones(N,1) / sqrt(N); % 根據(jù)虛警率計(jì)算門限 lambda Qinv(Pfa); %% 蒙特卡洛仿真 for k 1:length(SNR_dB) snr 10^(SNR_dB(k)/10); amp sqrt(snr); % 信號(hào)幅度保證信號(hào)能量 SNR T0 zeros(MC,1); % H0 下統(tǒng)計(jì)量 T1 zeros(MC,1); % H1 下統(tǒng)計(jì)量 for m 1:MC w randn(N,1); % 標(biāo)準(zhǔn)高斯白噪聲 x0 w; % H0純?cè)肼?x1 amp * s w; % H1信號(hào)加噪聲 T0(m) x0(:) * s; % 投影到信號(hào)方向 T1(m) x1(:) * s; end Pfa_sim(k) mean(T0 lambda); Pd_sim(k) mean(T1 lambda); end邏輯說明外層循環(huán)遍歷 SNR內(nèi)層循環(huán)做蒙特卡洛。amp * s構(gòu)造信號(hào)時(shí)由于s是單位能量向量amp的平方就是 SNR。門限lambda只依賴Pfa不隨 SNR 變化這正體現(xiàn)了參數(shù)已知 GLRT 的恒虛警特性。mean(T0 lambda)是統(tǒng)計(jì)超過門限的樣本比例用來驗(yàn)證虛警率是否落在 Pfa 附近mean(T1 lambda)就是檢測(cè)概率。參數(shù)說明MC取 100000 時(shí)虛警率估計(jì)的方差約為sqrt(Pfa*(1-Pfa)/MC)對(duì)Pfa1e-3約為 0.0001足夠看出和理論值的偏差。如果減小MC曲線會(huì)抖動(dòng)建議不低于 20000。仿真結(jié)束后不要急著畫圖先打印Pfa_sim和Pfa的偏差。如果偏差超過 20%多半是門限計(jì)算用了norminv(1-Pfa)而沒有注意左右尾方向或者統(tǒng)計(jì)量沒有做能量歸一化。4. 仿真結(jié)果解讀與性能評(píng)估4.1 結(jié)果.bmp 中關(guān)鍵曲線的判讀結(jié)果.bmp是這份源碼運(yùn)行后的輸出圖通常內(nèi)容因不同版本而異但最常見的形式是“檢測(cè)概率 Pd 隨 SNR 變化”的曲線圖橫軸 SNR_dB縱軸檢測(cè)概率同時(shí)可能疊加一條理論曲線。理論檢測(cè)概率的計(jì)算方式如下Pd_theory Q(lambda - sqrt(10.^(SNR_dB/10)));這條公式來自 H1 統(tǒng)計(jì)量 T ~ N(√SNR, 1) 超過門限 λ 的概率。對(duì)比仿真曲線和理論曲線重點(diǎn)看三點(diǎn)低 SNR 區(qū)域Pd 小于 0.2是否貼合這反映門限計(jì)算是否正確中間區(qū)域Pd 在 0.3~0.9是否有水平偏移這反映信號(hào)構(gòu)造的能量歸一化是否出錯(cuò)高 SNR 區(qū)域是否平滑收斂到 1這里蒙特卡洛樣本量不足會(huì)導(dǎo)致尾部抖動(dòng)。如果結(jié)果.bmp里還有第二條曲線比如經(jīng)驗(yàn)虛警率隨 SNR 變化的平坦線那說明仿真同時(shí)驗(yàn)證了恒虛警特性。理想情況下 Pfa_sim 應(yīng)該是一條水平直線幅度在 Pfa 附近隨機(jī)波動(dòng)。一旦發(fā)現(xiàn) Pfa_sim 隨 SNR 明顯傾斜就要回頭檢查噪聲產(chǎn)生是否用了有色噪聲或者信號(hào)歸一化時(shí)混入了噪聲功率。4.2 檢測(cè)概率與虛警概率的權(quán)衡參數(shù)已知 GLRT 沒有自由調(diào)節(jié)的“靈敏度旋鈕”唯一能改的就是虛警率 Pfa。虛警率每降一個(gè)數(shù)量級(jí)門限大約右移 0.5~1導(dǎo)致檢測(cè)概率曲線整體右移。下表給出了在同一 SNR 下不同 Pfa 對(duì) Pd 的影響數(shù)據(jù)基于 N10、SNR6 dB 的理論計(jì)算Pfa門限 λPd Q(λ - sqrt(SNR))1e-22.32630.9051e-33.09020.8401e-43.71900.7681e-54.26490.691可以看到為了把虛警從 1e-2 壓到 1e-5檢測(cè)概率損失約 0.2。工程上選擇 Pfa 時(shí)要考慮后續(xù)處理的代價(jià)雷達(dá)一次掃描虛警太多會(huì)引起計(jì)算機(jī)飽和通信系統(tǒng)誤同步會(huì)導(dǎo)致整個(gè)數(shù)據(jù)包報(bào)廢。我一般習(xí)慣先定虛警率上界再反推所需 SNR用這個(gè)仿真代碼先跑一遍理論曲線確認(rèn)系統(tǒng)預(yù)算夠不夠而不是直接調(diào)門限碰運(yùn)氣。4.3 常見調(diào)試問題與坑第一個(gè)坑是統(tǒng)計(jì)量方向?qū)懛础 * s和s * x在實(shí)數(shù)域結(jié)果相同但如果信號(hào)里含有復(fù)數(shù)分量忘記取共軛轉(zhuǎn)置就會(huì)得到錯(cuò)誤的投影。建議統(tǒng)一寫成x(:) * s(:)并且用(x * s)的實(shí)部作為判決量因?yàn)樘摬恐回暙I(xiàn)噪聲。第二個(gè)坑是門限單位不一致。Q 函數(shù)逆給出的分位數(shù)是標(biāo)準(zhǔn)正態(tài)尺度而有些實(shí)現(xiàn)里門限直接設(shè)為amp^2/2之類的能量閾值兩者差了sqrt(SNR)倍。判斷方法是把 SNR 設(shè)為 0 dB如果檢測(cè)概率接近 0.5 而虛警率不等于 Pfa說明統(tǒng)計(jì)量或門限沒有對(duì)準(zhǔn)尺度。第三個(gè)坑和 Qinv 的數(shù)值特性有關(guān)。當(dāng)Pfa小于 1e-8 時(shí)erfcinv的參數(shù)無(wú)限接近 2雙精度計(jì)算會(huì)丟失有效數(shù)字。此時(shí)我一般改用對(duì)數(shù)域近似% 極小虛警率的對(duì)數(shù)域近似適用于 Pfa 1e-8 t sqrt(-2 * log(Pfa) - log(2*pi) - 2*log(lambda));這個(gè)近似在 Pfa 很小時(shí)誤差可以忽略但不建議在常規(guī)范圍內(nèi)使用因?yàn)榈蠼飧菀滓肴藶檎`差。仿真時(shí)如果發(fā)現(xiàn) Pfa_sim 和設(shè)定值系統(tǒng)性偏差優(yōu)先檢查是否觸發(fā)了這種數(shù)值邊界。5. 把 GLRT 仿真擴(kuò)展到實(shí)際應(yīng)用場(chǎng)景5.1 參數(shù)失配時(shí)的魯棒性驗(yàn)證實(shí)際工程中“參數(shù)已知”往往是理想假設(shè)。為了評(píng)估失配的影響可以改造主流程里的信號(hào)向量比如發(fā)射端信號(hào)是s_true接收端匹配濾波器用的是s_hyp兩者之間存在頻率偏移或多普勒失配。做法是在仿真循環(huán)里用s_true生成數(shù)據(jù)用s_hyp計(jì)算統(tǒng)計(jì)量s_true exp(1j*2*pi*fd*(0:N-1)); % 真實(shí)信號(hào)帶多普勒頻移 s_hyp ones(N,1)/sqrt(N); % 假設(shè)信號(hào)無(wú)頻偏 % 生成數(shù)據(jù)用 s_true計(jì)算統(tǒng)計(jì)量用 s_hyp T abs(x * s_hyp); % 幅度檢測(cè)相位未知時(shí)取模參數(shù)說明fd是多普勒頻移與采樣間隔的乘積。失配會(huì)導(dǎo)致輸出 SNR 下降具體損耗因子是abs(s_true * s_hyp)^2 / (N^2)。你可以畫一條“檢測(cè)概率 vs 頻偏”的曲線觀察 GLRT 的性能懸崖在哪里。這個(gè)實(shí)驗(yàn)比單純跑理想曲線更有工程價(jià)值因?yàn)橥秸`差是真實(shí)系統(tǒng)必然存在的。5.2 用蒙特卡洛仿真替換固定閾值標(biāo)準(zhǔn) GLRT 的閾值由高斯假設(shè)解析給出但實(shí)際噪聲可能帶重尾或存在脈沖干擾。此時(shí)解析門限不再可靠更穩(wěn)妥的做法是先在純?cè)肼晽l件下跑大量蒙特卡洛得到經(jīng)驗(yàn)分布的 1-Pfa 分位數(shù)用這個(gè)分位數(shù)作為門限。實(shí)現(xiàn)如下MC_thresh 1000000; T_noise zeros(MC_thresh,1); for m 1:MC_thresh w randn(N,1) 0.2 * randn(N,1).^3; % 示例重尾噪聲 T_noise(m) w * s; end lambda_emp quantile(T_noise, 1 - Pfa);這段代碼的要點(diǎn)在于quantile得到的是經(jīng)驗(yàn)門限不需要假設(shè)噪聲閉式分布。代價(jià)是計(jì)算量增大但一旦噪聲模型變化只需重新跑一次T_noise而檢測(cè)部分的蒙特卡洛可以沿用同一個(gè)門限。注意MC_thresh要比1/Pfa至少大一個(gè)數(shù)量級(jí)否則尾部分位數(shù)估計(jì)偏差很大。5.3 從檢測(cè)到估計(jì)結(jié)合 MLE 的進(jìn)階思路參數(shù)已知 GLRT 的下一步自然擴(kuò)展是把信號(hào)幅度或相位作為未知量處理。做法是把檢測(cè)統(tǒng)計(jì)量里的s替換成最大似然估計(jì)得到的s_hat。以幅度未知、波形已知為例幅度 MLE 是a_hat (x*s)/(s*s)對(duì)應(yīng)的 GLRT 統(tǒng)計(jì)量變?yōu)門_glrt abs(x*s)^2 / (s*s); % 能量型統(tǒng)計(jì)量這個(gè)量在 H0 下服從指數(shù)分布H1 下服從非中心卡方分布。你可以在現(xiàn)有代碼里加一個(gè)分支比較參數(shù)已知和幅度未知兩種情況下檢測(cè)概率的差異通常幅度未知會(huì)帶來約 1.5~2 dB 的損失。這個(gè)擴(kuò)展不需要改動(dòng)Qinv.m只需把門限換成卡方分布的分位數(shù)chi2inv(1-Pfa, 1)。這樣一套仿真代碼就同時(shí)覆蓋了參數(shù)已知與參數(shù)部分未知兩類問題后續(xù)做自適應(yīng)檢測(cè)時(shí)可以直接復(fù)用這里的噪聲生成和蒙特卡洛框架。本文還有配套的精品資源點(diǎn)擊獲取