雷达目标检测之恒虚警检测（CA/SOCA/OS/GOCA-CFAR）【附MATLAB代码】

文章来源：微信公众号：EW Frontier

MATLAB主代码:

clear; close all; clc;
% 生成一个例子，然后用以下算法处理
% CA CFAR, OS CFAR, SOCA CFAR, GOCA CFAR, S-CFAR, Log CFAR
% 位置在50和55处有目标，50处信噪比10dB，55处信噪比15dB
% 前100单元噪声功率20dB，后100单元噪声功率30dB
 
num_cell = 200;
Pfa = 10^(-5);
 
signal1 = generateDataGaussianWhite(100, [50, 55], [35, 40], 20);
signal2 = generateDataGaussianWhite(100, [], [], 30);
signal = [signal1, signal2];
plot(1:num_cell, pow2db(signal), 'k-', 'linewidth', 0.5);
hold on;
 
% CA CFAR处理
[position, threshold, start_cell, stop_cell] = cacfar(signal, Pfa, 10, 2);
plot(start_cell:stop_cell, pow2db(threshold), 'linewidth', 1);
% OS CFAR处理
[position, threshold, start_cell, stop_cell] = oscfar(signal, Pfa, 10, 2, 15);
plot(start_cell:stop_cell, pow2db(threshold), 'linewidth', 1);
% SOCA CFAR处理
[position, threshold, start_cell, stop_cell] = socacfar(signal, Pfa, 10, 2);
plot(start_cell:stop_cell, pow2db(threshold), 'linewidth', 1);
% GOCA CFAR处理
[position, threshold, start_cell, stop_cell] = gocacfar(signal, Pfa, 10, 2);
plot(start_cell:stop_cell, pow2db(threshold), 'linewidth', 1);
 
 
grid on;
legend('信号', 'CA CFAR阈值', 'OS CFAR阈值', 'SOCA CFAR阈值', 'GOCA CFAR阈值');

CA-CFAR

function [position, threshold, start_cell, stop_cell] = cacfar(signal, Pfa, ref_num, guard_num)
    % ======>INPUT:
    % signal: Data of signal(include signal and noise).[DATATYPE: row vector]
    % Pfa: Probability of false alarm.[DATATYPE: scalar]
    % ref_num: Number of reference cell.[DATATYPE: scalar]
    % guard_num: Number of guard cell.[DATATYPE: scalar]
    % ======>OUTPUT:
    % position: positions of target.[DATATYPE: row vector]
    % threshold: CFAR threshold of input signal.[DATATYPE: row vector]
    position = [];
    left_num = guard_num + ref_num;
    start_cell = left_num + 1;
    stop_cell = length(signal) - left_num;
    N = 2*ref_num;
    alpha = N * (Pfa ^ (-1/N) - 1);
 
    threshold = zeros(1, stop_cell - start_cell + 1);
    for ii = start_cell : stop_cell
        tmp_data = [signal(ii-left_num : ii-guard_num-1), ...
                    signal(ii+guard_num+1 : ii+left_num)];
        tmp = mean(tmp_data) * alpha;
        % threshold(ii - left_num)
        % tmp
        threshold(ii - left_num) = tmp;
        if tmp < signal(ii)
            position = [position, ii];
        end
    end
return

OS-CFAR

function [position, threshold, start_cell, stop_cell] = oscfar(signal, Pfa, ref_num, guard_num, k)
    % ======>INPUT:
    % signal: Data of signal(include signal and noise).[DATATYPE: row vector]
    % Pfa: Probability of false alarm.[DATATYPE: scalar]
    % ref_num: Number of reference cell.[DATATYPE: scalar]
    % guard_num: Number of guard cell.[DATATYPE: scalar]
    % k: # k sorted statistical value.[DATATYPE: scalar]
    % ======>OUTPUT:
    % position: positions of target.[DATATYPE: row vector]
    % threshold: CFAR threshold of input signal.[DATATYPE: row vector]
    position = [];
    left_num = guard_num + ref_num;
    start_cell = left_num + 1;
    stop_cell = length(signal) - left_num;
    N = 2*ref_num;
    alpha = get_alpha(Pfa, N, k);
 
    threshold = zeros(1, stop_cell - start_cell + 1);
    for ii = start_cell : stop_cell
        tmp_data = [signal(ii-left_num : ii-guard_num-1), ...
                    signal(ii+guard_num+1 : ii+left_num)];
        sorted_data = sort(tmp_data);
        T = sorted_data(k) * alpha;
        threshold(ii - left_num) = T;
        if T < signal(ii)
            position = [position, ii];
        end
    end
 
    function alpha = get_alpha(Pfa_set, N, k)
        % 用中点法求alpha_OS
        left_alpha = 0;
        while true
            right_alpha = left_alpha + 1;
            this_pfa = k * nchoosek(N, k) * beta(right_alpha + N - k + 1, k);
            if this_pfa < Pfa_set
                break;
            end
            left_alpha = right_alpha;
        end
        
        mid_alpha = 0.5 * (left_alpha + right_alpha);
        this_pfa = k * nchoosek(N, k) * beta(mid_alpha + N - k + 1, k);
        while abs(this_pfa - Pfa_set) > 0.000001 * Pfa_set
            if this_pfa > Pfa_set
                left_alpha = mid_alpha;
            else
                right_alpha = mid_alpha;
            end
            mid_alpha = 0.5 * (left_alpha + right_alpha);
            this_pfa = k * nchoosek(N, k) * beta(mid_alpha + N - k + 1, k);
        end
        % this_pfa
        alpha = mid_alpha;
    return
return

MATLAB仿真结果：