GLRT_AMF_ACE程序源码matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MF的阈值不知 % 名称：Kelly's GLRT and AMF and ACE and Rao % 时间：2017.09.13 % 作者：王莎 % 功能：仿真检测概率Pd与信噪比SNR的关系曲线 % 注意：变量命名中1代表GLRT,2代表AMF,3代表ACE,4代表Rao % AMF ACE的检测门限需要利用划分门限逼近预设虚警率的方法得到 % 当图不正常时，首先想到是不是计算AMF ACE给定的阈值范围太小 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc clear all close all %% 常数定义 deg2rad = pi / 180; % deg -> rad rad2deg = 180 / pi; % rad -> deg ima = sqrt(-1); % 定义虚数单位 % 仿真参数设置(修改时必须同时修改ComputeAMFpfa.m和ComputeACEpfa.m中的参数) Na = 16; % 阵元数 Np = 1; % 脉冲串数 N = Na * Np; % 主采样数据维数 K = 2 * N; % 辅助采样数据个数 PFA = 1.0e-6; % 预设虚警率 L = K + 1 - N; % 阵列参数设置 f = 1.5e3; % 单频信号频率（Hz） c = 1500; % 介质中声速（m/s）ou lamda = c * (1/f); % 声波波长(m) dd = 1/2 * lamda; % 线阵阵元间距(m) d=0 : dd : (N-1)*dd; % 每个阵元所在的位置 ang_deg = 0; % 期望信号角(角度) ang_rad = ang_deg * deg2rad; % 期望信号角(弧度) omiga_t = 0.4; % 目标Doppler %目标空时信号 aN_t = zeros(Na,1); bN_t = zeros(Np,1); aN_t = exp(-ima * pi * [0:Na-1]' * cos(ang_rad)) / sqrt(Na); % 空域导向向量 bN_t = exp(-ima * pi * [0:Np-1]' * omiga_t) / sqrt(Np); % 时域导向向量 ss = zeros(N,1); %空时导向向量 ss = kron(aN_t,bN_t); %% 计算检测门限 % GLRT(利用公式) Ratio1 = 1 / PFA^(1/L); % 似然比值，利用公式PFA = (1/Ratio)^(K+1-N) Threshold1 = (Ratio1 - 1) / Ratio1;% 似然比检测门限 %AMF 计算虚警率所对应的门限(带入虚警率与门限满足的等式中) ThresholdRange2 = 0:0.0001:10; % 用于确定门限的位置，最大范围(10)有时需要扩大，以0.0001为单位保证精度 len2 = length(ThresholdRange2); pfa_real2 = zeros(1,len2); for i = 1:length(ThresholdRange2) pfa_real2(1,i) = ComputeAMFpfa(ThresholdRange2(i),N,K); if pfa_real2(1,i)-PFA<eps Threshold2 = ThresholdRange2(i); break; end end % 另一种方法：解非线性方程，但与初值有很大关系。不太正确 % a0 = 0.5 ;%阈值初值迭代，GLRT的初值 % Threshold = fsolve(@myequs,a0,optimset('Display','on'));%检测门限 % ACE 计算虚警率所对应的门限(带入虚警率与门限满足的等式中) ThresholdRange3 = 0:0.001:20; % 用于确定门限的位置 len3 = length(ThresholdRange3); pfa_real3 = zeros(1,len3); Threshold = 0; for i = 1:length(ThresholdRange3) pfa_real3(1,i) = ComputeACEpfa(ThresholdRange3(i),N,K); if pfa_real3(1,i)-PFA<eps Threshold = ThresholdRange3(i); break; end end Threshold3 = Threshold/(1+Threshold); % Rao(利用公式) Threshold4 = 1 - PFA^(1/K); %% 构造s,n与zk % X = diag(rand(N,1)); % U = orth(rand(N,N)); % M = U*X*U'; % 理想协方差矩阵(正定对称矩阵) % [Ev,D,En] = svd(M); c = linspace(10,1,N); % 混响的协方差矩阵为Toeplitz矩阵，实际中多为混响忽略噪声 r = linspace(10,1,N); % 用于产生Toeplitz矩阵 M = toeplitz(c,r); % 构造理想协方差矩阵 [Ev,D,En] = svd(M); % M = UDV'得到信号子空间Ev和噪声子空间En,对角阵D SNR = 5 : 0.5 : 25; % 信噪比取值范围 T = length(SNR); Pd1 = zeros(T,1); % GLRT检测概率矩阵 Pd2 = zeros(T,1); % AMF检测概率矩阵 Pd3 = zeros(T,1); % ACE检测概率矩阵 Pd4 = zeros(T,1); % Rao检测概率矩阵 s = M^(-1/2) * ss; % 白化处理 for j = 1 : length(SNR) num1 = 0; % GLRT同一SNR下超过阈值的次数 num2 = 0; % AMF同一SNR下超过阈值的次数 num3 = 0; % ACE同一SNR下超过阈值的次数 num4 = 0; % Rao同一SNR下超过阈值的次数 Count = 1e4; % 蒙特卡洛仿真次数 for i = 1 : Count n = M^(-1/2) * En * (randn(N,1) + ima * randn(N,1)) / sqrt(2); % 主采样z = bs+n 中的噪声n zk = M^(-1/2) * En * (randn(N,K) + ima * randn(N,K)) / sqrt(2); % 由噪声子空间生成辅助数据，均值为0的高斯噪声 S = zk * zk'; % NoiseCov = S / K; % 由K个辅助数据来估计噪声协方差矩阵,zk'表示共轭转置 % NoisePwr = s' / NoiseCov * s; % SNR = abs(b)^2*(s'*NoiseCov*s) b = sqrt(10^(SNR(j)/10)); z = b * s + n; % 主采样数据 y1 = abs(s'/S*z)^2 / ((s'/S*s) * (1+(z'/S*z))); % GLRT检测统计量 y2 = abs(s'/S*z)^2 / (s'/S*s); % AMF检测统计量 y3 = y2 / (z'/S*z); % ACE检测统计量 y4 = (y1^2) / (y2*(1-y1)); % Rao检测统计量 if y1>Threshold1 num1 = num1 + 1; end if y2>Threshold2 num2 = num2 + 1; end if y3>Threshold3 num3 = num3 + 1; end if y4>Threshold4 num4 = num4 + 1; end end Pd1(j,1) = num1 / Count; % GLRT同一SNR下的Pd Pd2(j,1) = num2 / Count; % AMF同一SNR下的Pd Pd3(j,1) = num3 / Count; % ACE同一SNR下的Pd Pd4(j,1) = num4 / Count; % Rao同一SNR下的Pd end % figure % plot(SNR,Pd1,'-b',SNR,Pd2,'-r') % legend('GLRT','AMF') % xlabel('SNR(/dB)'),ylabel('Pd') % text(21,0.15,'N=10,K=20') % text(21,0.08,'PFA=1.0 E-6') % title('GLRT AMF检测概率与信噪比的关系') % figure % plot(SNR,Pd1,'-b',SNR,Pd2,'-r',SNR,Pd3,'-k') % legend('GLRT','AMF','ACE') % xlabel('SNR(/dB)'),ylabel('Pd') % text(21,0.15,'N=10,K=20') % text(21,0.08,'PFA=1.0 E-6') % title('GLRT AMF ACE检测概率与信噪比的关系') figure plot(SNR,Pd1,'-b',SNR,Pd2,'-r',SNR,Pd3,'-k',SNR,Pd4,'-g') legend('GLRT','AMF','ACE','Rao') xlabel('SNR(/dB)'),ylabel('Pd') text(21,0.15,'N=16,K=32') text(21,0.08,'PFA=1.0 E-6') title('GLRT AMF ACE Rao检测概率与信噪比的关系')