www.pudn.com > radarsignalprocessing.zip > radar_signal_detection.m, change:2013-12-11,size:12240b


%% 雷达信号检测算法仿真实例,包括MTI、MTD、CFAR等处理。 
 
clear all; 
close all; 
clc; 
 
%% 动目标 
 
fr = 1e3;                                                                  % 脉冲重复频率(Hz) 
tr = 1 / fr;                                                               % 脉冲重复间隔(s) 
fs = 1e6;                                                                  % 快时间维采样频率(Hz)(等于发射脉冲带宽,发射脉冲带宽通常为载波频率的1%) 
ts = 1 / fs;                                                               % 快时间维采样间隔(s) 
fc = 100e6;                                                                % 载波频率(Hz) 
vr1 = 100;                                                                 % 动目标1相对雷达的径向速度(m/s)   
vr2 = 700;                                                                 % 动目标2相对雷达的径向速度(m/s)   
c = 3e8;                                                                   % 电磁波速(m/s)                                          
lamada = c / fc;                                                           % 载波波长(m) 
fd1 = 2 * vr1 / lamada;                                                    % 动目标1的多普勒频率(Hz)(或多普勒移动,即发射频率与接收频率之差) 
fd2 = 2 * vr2 / lamada;                                                    % 动目标2的多普勒频率(Hz)(或多普勒移动,即发射频率与接收频率之差) 
 
M = fs / fr;                                                               % 1个脉冲重复间隔内的采样点数(快时间维的距离点数) 
M10 = 10 * M;                                                              % 10个脉冲重复间隔内的采样点数 
m = 1 : M10; 
target1 = 20 * exp(j * (2 * m * pi * fd1 / fr));                           % 动目标1 
target2 = 20 * exp(j * (2 * m * pi * fd2 / fr));                           % 动目标2 
 
figure(1); 
subplot(211); 
plot(real(target1)); 
title('动目标1——实部'); 
xlabel('距离点'); 
grid on; 
subplot(212); 
plot(imag(target1)); 
title('动目标1——虚部'); 
xlabel('距离点'); 
grid on; 
figure(2); 
subplot(211); 
plot(real(target2)); 
title('动目标2——实部'); 
xlabel('距离点'); 
grid on; 
subplot(212); 
plot(imag(target2)); 
title('动目标2——虚部'); 
xlabel('距离点'); 
grid on; 
f1 = linspace(0, fr, length(target1)); 
figure(3); 
plot(f1, abs(fft(target1))); 
title('动目标1幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('幅度'); 
grid on; 
f2 = linspace(0, fr, length(target2)); 
figure(4); 
plot(f2, abs(fft(target2))); 
title('动目标1幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('幅度'); 
grid on; 
 
%% 杂波 
 
randn('state', sum(100 * clock)); 
clutteri = 5 * randn(1, length(target1)); 
randn('state', sum(10 * clock)); 
clutterq = 5 * randn(1, length(target1)); 
clutter = clutteri + j * clutterq;                                         % 杂波 
 
figure(5); 
subplot(211); 
plot(real(clutter)); 
title('杂波——实部'); 
xlabel('距离点'); 
grid on; 
subplot(212); 
plot(imag(clutter)); 
title('杂波——虚部'); 
xlabel('距离点'); 
grid on; 
f3 = linspace(0, fr, length(clutter)); 
figure(6); 
plot(f2, abs(fft(clutter))); 
title('杂波幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('幅度'); 
grid on; 
 
%% 动目标与杂波的混合信号 
 
s = target1 + target2 + clutter;                                           % 动目标与杂波的混合信号 
 
figure(7); 
subplot(211); 
plot(real(s)); 
title('动目标与杂波的混合信号——实部'); 
xlabel('距离点'); 
grid on; 
subplot(212); 
plot(imag(s)); 
title('动目标与杂波的混合信号——虚部'); 
xlabel('距离点'); 
grid on; 
f5 = linspace(0, fr, length(s)); 
figure(8); 
plot(f5, abs(fft(s))); 
title('混合信号幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('幅度'); 
grid on; 
 
%% MTI杂波抑制(三脉冲对消处理) 
 
%————————————————对消运算————————————————— 
x = ones(1, 4);                                                            % 动目标所在的距离点 
n = length(x); 
points1 = [zeros(1, 300), x, zeros(1, 600), x, zeros(1, M - 900 - 2 * n)]; % 1个脉冲重复间隔内M个距离点(将动目标1所在的距离点301、302、303、304;动目标2所在的距离点905、906、907、908与杂波所在距离点分开) 
points10 = repmat(points1, 1, 10);                                         % 10个脉冲重复间隔内10*M个距离点 
 
s2 = points10 .* (target1 + target2) + clutter;                            % 动目标与杂波的混合信号(确定动目标所在的距离点301、302、303、304、905、906、907、908,并为单元平均CFAR检测的假设提供前提,即动目标单纯占据8个距离点,其余距离点为杂波) 
 
mti_in = reshape(s2, length(s2) / 10, 10)';                                % 形成由10个脉冲组成的二维数据矩阵(慢时间维10个脉冲周期,快时间维1000个距离点) 
for t = 1 : 8 
    mti_out(t, :) = mti_in(t, :) - 2 * mti_in(t + 1, :) + mti_in(t + 2, :);% 三脉冲对消运算(同一距离点在3个相邻脉冲重复间隔内的数据做相减运算) 
end 
 
%———————————二脉冲与三脉冲对消器的频率响应——————————— 
f5 = -3 * fr : 0.01 : 3 * fr; 
H2 = 1 - exp(-1 * j * 2 * pi * f5 / fr);                                   % 二脉冲对消器的频率响应 
H3 = 1 - 2 * exp(-1 * j * 2 * pi * f5 / fr) + exp(-2 * j * 2 * pi * f5 / fr);  % 三脉冲对消器的频率响应 
 
figure(9); 
plot(f5, abs(H2), 'b', f5, abs(H3), 'r'); 
title('二脉冲和三脉冲对消器的幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('频响幅度'); 
legend('二脉冲', '三脉冲'); 
grid on; 
figure(10); 
plot(f5, 20 * log10(abs(H2)), 'b', f5, 20 * log10(abs(H3)), 'r'); 
title('二脉冲和三脉冲对消器的幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('频响幅度/dB'); 
legend('二脉冲', '三脉冲'); 
grid on; 
 
%————————————显示三脉冲对消杂波抑制效果———————————— 
b = [1, -2, 1];                                                            % 三脉冲对消器系统函数分子多项式系数向量 
a = [1];                                                                   % 三脉冲对消器系统函数分母多项式系数向量 
y = filter(b, a, s);                                                       % 混合信号s经过三脉冲对消器输出信号y 
 
f5 = linspace(0, fr, length(y)); 
figure(11); 
plot(f5, abs(fft(y))); 
title('三脉冲对消器输出信号幅频特性'); 
xlabel('频率(Hz)'); 
ylabel('幅度'); 
grid on; 
 
%% 脉冲多普勒处理(多普勒谱分析,也可理解为MTD多通道滤波处理) 
 
K = 8;                                                                     % FFT采样点数即MTD通道数 
for m = 1 : M 
    mtd_out(:, m) = fft(mti_out(:, m) .* hamming(K), K);                   % 对同一距离点不同脉冲重复间隔的数据做K点FFT运算(含加窗) 
end 
p = 1 : 1 : M;                                                             % 距离点 
 
figure(12); 
plot(p, abs(mtd_out(1, :))); 
title('脉冲多普勒处理第0通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(13); 
plot(p, abs(mtd_out(2, :))); 
title('脉冲多普勒处理第1通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(14); 
plot(p, abs(mtd_out(3, :))); 
title('脉冲多普勒处理第2通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(15); 
plot(p, abs(mtd_out(4, :))); 
title('脉冲多普勒处理第3通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(16); 
plot(p, abs(mtd_out(5, :))); 
title('脉冲多普勒处理第4通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(17); 
plot(p, abs(mtd_out(6, :))); 
title('脉冲多普勒处理第5通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(18); 
plot(p, abs(mtd_out(7, :))); 
title('脉冲多普勒处理第6通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
figure(19); 
plot(p, abs(mtd_out(8, :))); 
title('脉冲多普勒处理第7通道结果'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on; 
 
%% 自适应门限检测(单元平均CFAR处理) 
 
nr = 8;                                                                    % 待检测距离点每侧的参考点数 
np = 3;                                                                    % 待检测距离点每侧的保护点数 
Pfa = 10e-3;                                                               % 虚警概率 
alafa = 2 * nr * (Pfa ^ (-1 / (2 * nr)) - 1);                              % 门限因子 
temp = zeros(K, M); 
 
for u = 1 : 8                                                              % 分别对每个多普勒采样点做检测判决 
    z = mtd_out(u, :); 
     
    for i = 1 : 1 + np                                                     % 第i个待检测距离点 
        sum = 0; 
        for t = i + 1 + np : i + 1 + nr                                    % 第t个参考距离点 
            sum = sum + (abs(z(t))) .^ 2;                                   
        end 
        itfPower(i) = sum / nr;                                            % 估计第i个待检测距离点的干扰功率 
    end 
     
    for i = 1 + 1 + np : np + nr 
        sum1 = 0; 
        for t = i * 2 - 1 : (i * 2 - 1) + nr - 1 
            sum1 = sum1 + (abs(z(t))) .^ 2;                                % 第i个待检测距离点的右侧参考距离点干扰功率估计 
        end 
        sum2 = 0; 
        for t = 1 : i - 1 - np 
            sum2 = sum2 + (abs(z(t))) .^ 2;                                % 第i个待检测距离点的左侧参考距离点干扰功率估计 
        end 
        itfPower(i) = (sum1 + sum2) / 2; 
    end 
     
    for i = 1 + np + nr : M - np - nr 
        sum4 = 0; 
        for t = i + 1 + np : i + 1 + nr; 
            sum4 = sum4 + (abs(z(t))) .^ 2; 
        end 
        sum3 = 0; 
        for t = i - np - nr : i - 1 - np 
            sum3 = sum3 + (abs(z(t))) .^ 2; 
        end 
        itfPower(i) = (sum4 + sum3) / 2; 
    end 
     
    for i = M - nr - np + 1 : M - np - 1 
        sum5 = 0; 
        for t = i + 1 + np : M 
            sum5 = sum5 + (abs(z(t))) .^ 2; 
        end 
        sum6 = 0; 
        for t = i - np - nr : i - 1 - np 
            sum6 = sum6 + (abs(z(t))) .^ 2; 
        end 
        itfPower(i) = (sum5 + sum6) / 2; 
    end 
     
    for i = M - np : M 
        sum7 = 0; 
        for t = i - nr - np : i - 1 - nr; 
            sum7 = sum7 + (abs(z(t))) .^ 2; 
        end 
        itfPower(i) = sum7 / nr; 
    end 
     
    s_cfar = zeros(1, M); 
     
    for i = 1 : M                                                           
        T(i) = alafa * itfPower(i);                                        % 第i个距离点的检测门限 
        if (abs(z(i))) .^ 2 > T(i)                                         % 判断第i个距离点的信号功率是否高于其检测门限 
            s_cfar(i) = z(i);                                              % 若高于检测门限,判决该距离点存在动目标,将其输出 
        end 
    end 
     
    temp(u, :) = s_cfar;                                                   % 将每个多普勒采样点的判决结果存储在临时矩阵中 
     
    figure(20); 
    subplot(2, 4, u); 
    plot(p, 10 * log10((abs(z)) .^ 2), 'r', p, 10 * log10(T), 'g'); 
    legend('混合信号', '检测门限', 4); 
    xlabel('距离点'); 
    ylabel('功率/dB'); 
    grid on; 
    figure(21); 
    subplot(2, 4, u); 
    plot(p, abs(s_cfar)); 
    xlabel('距离点'); 
    ylabel('幅度'); 
    grid on; 
end 
 
figure(20); 
subplot(241); 
title('CA CFAR 第0通道混合信号与检测门限'); 
subplot(242); 
title('CA CFAR 第1通道混合信号与检测门限'); 
subplot(243); 
title('CA CFAR 第2通道混合信号与检测门限'); 
subplot(244); 
title('CA CFAR 第3通道混合信号与检测门限'); 
subplot(245); 
title('CA CFAR 第4通道混合信号与检测门限'); 
subplot(246); 
title('CA CFAR 第5通道混合信号与检测门限'); 
subplot(247); 
title('CA CFAR 第6通道混合信号与检测门限'); 
subplot(248); 
title('CA CFAR 第7通道混合信号与检测门限'); 
figure(21); 
subplot(241); 
title('CA CFAR 第0通道检测判决结果'); 
subplot(242); 
title('CA CFAR 第1通道检测判决结果'); 
subplot(243); 
title('CA CFAR 第2通道检测判决结果'); 
subplot(244); 
title('CA CFAR 第3通道检测判决结果'); 
subplot(245); 
title('CA CFAR 第4通道检测判决结果'); 
subplot(246); 
title('CA CFAR 第5通道检测判决结果'); 
subplot(247); 
title('CA CFAR 第6通道检测判决结果'); 
subplot(248); 
title('CA CFAR 第7通道检测判决结果'); 
 
target_cfar = zeros(1, M); 
for i = 1 : M 
    target_cfar(:, i) = abs(temp(1, i)) + abs(temp(2, i)) + abs(temp(3, i)) + abs(temp(4, i)) + abs(temp(5, i)) + abs(temp(6, i)) + abs(temp(7, i)) + abs(temp(8, i));                                    
                                                                           % 每个距离点取8个通道之合作为动目标信号的最终判决结果 
end 
 
figure(22); 
plot(p, target_cfar); 
title('CA CFAR 8通道合成检测结果——动目标信号'); 
xlabel('距离点'); 
ylabel('幅度'); 
grid on;