雷达无线电系列（三）经典CFAR算法门限因子alpha计算（matlab）

前言 本文汇集CA、SO、GO、OS、杂波图等恒虚警算法的门限因子求解方法及其函数

1，CA-CFAR 【非常简单，可以直接求解】

%% 均值恒虚警_门限因子计算公式
%% 版本:v1
%% 时间:2019.11.08
%% 终版【不在优化】

function [ alpha ] = form_ALPHA_ca( PFA,N )
%FORM_ALPHA_CA 此处显示有关此函数的摘要
%   PFA:虚警概率
%   N:参考单元个数

alpha=N.*(PFA.^(-1./N)-1);

end

　　

2，SO-CFAR　　【这个门限因子求解是一个非常复杂的高阶函数，基本上无法通过反函数求解，这里使用二分求解法解决】

PFA计算函数

%% 虚警概率_均值选小
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ PFA ] = form_PFA_so( ALPHA,N )
%FORM_PFA_GO 此处显示有关此函数的摘要
%   ALPHA:门限因子
%   N:参考窗个数

PFA=0;
n=N/2;
for i=0:n-1
    PFA=PFA+2.*gamma(n+i)./gamma(i+1)./gamma(n).*(2+ALPHA./n).^(-(n+i));
end

end

　　

2，GO-CFAR

PFA计算函数

%% 虚警概率_均值选大
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ PFA ] = form_PFA_go( ALPHA,N )
%FORM_PFA_GO 此处显示有关此函数的摘要
%   ALPHA:门限因子
%   N:参考窗个数

PFA=0;
n=N/2;
for i=0:n-1
    PFA=PFA-2.*gamma(n+i)./gamma(i+1)./gamma(n).*(2+ALPHA./n).^(-(n+i));
end
PFA=PFA+2.*(1+ALPHA./n).^(-n);

end

　　

3，OS-CFAR 

%% 虚警概率_有序
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ PFA ] = form_PFA_os( ALPHA,N,Rate )
%FORM_PFA_GO 此处显示有关此函数的摘要
%   ALPHA:门限因子
%   N:参考窗个数
%   Rate:比例点

k=ceil(N.*Rate);
PFA=gamma(N+1).*gamma(N-k+ALPHA+1)./gamma(N-k+1)./gamma(N+ALPHA+1);

end

　　

4，杂波图

点参数

%% 虚警概率_杂波图_点参数
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ PFA ] = form_PFA_cm_point( ALPHA,m,r )
%FORM_PFA_CM 此处显示有关此函数的摘要
%   ALPHA:门限因子
%   m:天线旋转周期

PFA=1;
for n=0:m-1
    PFA=PFA.*(1+ALPHA.*r.*(1-r).^n).^-1;
end

end

面参数

%% 虚警概率_杂波图_面参数
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ PFA ] = form_PFA_cm_surface( ALPHA,m,r,M )
%FORM_PFA_CM 此处显示有关此函数的摘要
%   ALPHA:门限因子
%   m:天线旋转周期
%   M:参考单元数

PFA=1;
for n=0:m-1
    PFA=PFA.*(1+ALPHA.*r.*(1-r).^n./M).^-1;
end
end

　　

5，二分求解法核心函数

%% 二分法求解方程的解
%% 时间：2019.11.27
%% 版本：v1【不在优化】

function [ d1 ] = func_SOLUTION_binary( d1_scope,d2_precision,func,parameter)
%FUNC_SOLUTION_BINARY 此处显示有关此函数的摘要
%   d1: 已知数值
%   d2: 目标数值
%   d2_scope: 目标初始范围【小值，大值】
%   d1_precision: 目标数值精度【精度值】
%   func: 公式
%   parameter: 相关参数【类型自定】 1=ALPHA, end=d2

shape=size(parameter);
if shape(2)==2      %只有一个参数额外参数需要输入
    while 1
        d1=mean(d1_scope);
        d2=func(d1,parameter(1,1));     % 本质上是PFA_pre
        d2_difference=1/d2-1/parameter(1,end);
        if abs(d2_difference)<d2_precision || abs(d1_scope(1,1)-d1_scope(1,2))<0.001
            return;
        elseif d2_difference<0
            d1_scope(1,1)=d1;
        else
            d1_scope(1,2)=d1;
        end
    end
elseif shape(2)==3      % 用于OS_CFAR门限因子的计算
    while 1
        d1=mean(d1_scope);
        d2=func(d1,parameter(1,1),parameter(1,2));     % 本质上是PFA_pre
        d2_difference=1/d2-1/parameter(1,end);
        if abs(d2_difference)<d2_precision || abs(d1_scope(1,1)-d1_scope(1,2))<0.001
            return;
        elseif d2_difference<0
            d1_scope(1,1)=d1;
        else
            d1_scope(1,2)=d1;
        end
    end
end
end

　　

6，测试样例及其输出结果

clear;

PFA=10^(-4);
N=36;
Rate=0.5;

[ ALPHA ] = form_ALPHA_ca( PFA,N );
[ PFA1 ] = form_PFA_ca( ALPHA,N );

d1_scope=[0,100];
d2_precision=1;
func1=@form_PFA_ca;
func2=@form_PFA_so;
func3=@form_PFA_go;
func4=@form_PFA_os;

func5=@form_PFA_cm_point;

parameter=[N,PFA];
parameter4=[N,Rate,PFA];
parameter5=[100,1./512,PFA];


[ ALPHA1 ] = func_SOLUTION_binary( d1_scope,d2_precision,func1,parameter);
[ ALPHA2 ] = func_SOLUTION_binary( d1_scope,d2_precision,func2,parameter);
[ ALPHA3 ] = func_SOLUTION_binary( d1_scope,d2_precision,func3,parameter);
ALPHA3_1=form_ALPHA_os (N,Rate,PFA);
[ ALPHA4 ] = func_SOLUTION_binary( d1_scope,d2_precision,func4,parameter4);

[ ALPHA5 ] = func_SOLUTION_binary( d1_scope,d2_precision,func5,parameter5);

输出结果展示：