文章来源:微信公众号: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仿真结果: