Matlab 恒虚警率检测CFAR

function [res, th_lr, th_ud]= cfar_filter(data,varargin)
% Function: 恒虚警率检测
%
% Input:
%   data: 数据
% Output:
%   res: 检测结果
%   th_lr:左右方向计算门限
%   th_ud:上下方向计算门限
%
% Example:
%     res = cfar(data);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

guard_r_lr = find_Para(varargin,'guard_r_lr',2); % 左右方向保护单元半径
training_r_lr = find_Para(varargin,'training_r_lr',10); % 左右方向训练单元半径
guard_r_ud = find_Para(varargin,'guard_r_ud',1); % 上下方向保护单元半径
training_r_ud = find_Para(varargin,'training_r_ud',2); % 上下方向训练单元半径
Pfa_lr = find_Para(varargin,'Pfa_lr',1e-10); % 左右方向虚警率
Pfa_ud = find_Para(varargin,'Pfa_ud',1e-4); % 上下方向虚警率
dim = find_Para(varargin,'dim',3); % 维度 1一维(左右方向) 2一维十字 3宽十字 4窗口
type = find_Para(varargin,'type',1); % 类型 1均值 2最小值 3最大值 4排序

row = size(data,1);
col = size(data,2);
res = false(size(data));
th_lr = zeros(size(data));
th_ud = zeros(size(data));

mid_row = training_r_ud+guard_r_ud+1;
mid_col = training_r_lr+guard_r_lr+1;
up_guard = training_r_ud + 1;
down_guard = training_r_ud + guard_r_ud * 2 + 1;
left_guard = training_r_lr + 1;
right_guard = training_r_lr + guard_r_lr * 2 + 1;

window_data_use_lr = false(training_r_ud*2+guard_r_ud*2+1,training_r_lr*2+guard_r_lr*2+1);
window_data_use_ud = false(training_r_ud*2+guard_r_ud*2+1,training_r_lr*2+guard_r_lr*2+1);
if dim >= 1
    window_data_use_lr(mid_row,:) = 1;
end
if dim >= 2
    window_data_use_ud(:,mid_col) = 1;
end
if dim >= 3
    window_data_use_lr(up_guard:down_guard,:) = 1;
    window_data_use_ud(:,left_guard:right_guard) = 1;
end
if dim >= 4
    window_data_use_lr(:,[1:left_guard-1,right_guard+1:end]) = 1;
end
window_data_use_lr(up_guard:down_guard,left_guard:right_guard) = 0;
window_data_use_ud(up_guard:down_guard,left_guard:right_guard) = 0;

Nc_lr = sum(sum(window_data_use_lr));
Nc_ud = sum(sum(window_data_use_ud));
CAThreshold_lr = Nc_lr*(Pfa_lr^(-1/Nc_lr)-1);
CAThreshold_ud = Nc_ud*(Pfa_ud^(-1/Nc_ud)-1);

data_frame = nan(row+(guard_r_ud+training_r_ud)*2,col+(guard_r_lr+training_r_lr)*2);
row_shift = guard_r_ud+training_r_ud;
col_shift = guard_r_lr+training_r_lr;
data_frame(1+row_shift:row+row_shift,1+col_shift:col+col_shift) = data;
for i = 1:row
    for j = 1:col
        window_data_lr = data_frame(i:i+row_shift*2,j:j+col_shift*2);
        window_data_ud = window_data_lr;
        window_data_lr(~window_data_use_lr) = nan;
        window_data_ud(~window_data_use_ud) = nan;

        th_lr(i,j) = get_th(window_data_lr,CAThreshold_lr,type);
        if dim == 1
            res(i,j) = data(i,j) > th_lr(i,j);
        else
            th_ud(i,j) = get_th(window_data_ud.',CAThreshold_ud,type);
            res(i,j) = data(i,j) > th_lr(i,j) & data(i,j) > th_ud(i,j);
        end
    end
end
end

function value = find_Para(var,name,value)
% var:varargin
% name: 参数名
% value: 默认值
if size(var,2) == 1
    var = var{:};
end
if ~isempty(var)
    ind = find(strcmp(var,name));
    if ind
        v = var{ind+1};
        try
            if ~isempty(v) && ~isnan(v)
                value = v;
            end
        catch
            value = v;
        end
    end
end
end

function th = get_th(data,CAThreshold,type)
data = data(~isnan(data));
len = round(length(data)/2);
if type == 1
    th = mean(data)*CAThreshold;
elseif type == 2
    th = min([mean(data(1:len)),mean(data(len+1:end))])*CAThreshold;
elseif type == 3
    th = max([mean(data(1:len)),mean(data(len+1:end))])*CAThreshold;
elseif type == 4
    sort_window_data = sort(data(:));
    th = sort_window_data(round(0.75 * length(sort_window_data)))*CAThreshold;
end
end
%% CFAR Test
% rs = RandStream('mt19937ar','Seed',2010);
% npower = db2pow(-10);  % Assume 10dB SNR ratio
% Npoints = 1e4;
% rsamp = randn(rs,Npoints,1)+1i*randn(rs,Npoints,1);
% ramp = linspace(1,10,Npoints)';
% x = abs(sqrt(npower*ramp./2).*rsamp).^2;
% [x_detected1,th] = cfar_filter(x','guard_r_lr',1,'training_r_lr',100,'Pfa_lr',1e-3,'dim',1,'type',1);
% act_pfa = sum(x_detected1)/Npoints
% figure
% plot(1:length(x),x,1:length(x),th,...
%     find(x_detected1),x(x_detected1),'o')
% legend('Signal','Threshold','Detections','Location','Northwest')
% xlabel('Time Index')
% ylabel('Level')