ssy的小天地

愿汝如是 千金不换

0%

FPGA数字信号处理之CFAR检测

发表于 更新于 分类于

本节主要介绍了一维CFAR的基本原理以及它的MATLAB仿真与Verilog硬件实现的思路。

一维CFAR的基本原理

  • 毫米波雷达中，CFAR（Constant False Alarm Rate）算法是一种常用的目标检测和跟踪算法。它的主要作用是在背景噪声（杂波）中检测出目标信号，同时保证误检概率不变。
  • CFAR算法的步骤如下：
    • 确定检测窗口的大小和形状，例如矩形、圆形等。
    • 在检测窗口内计算信号功率的平均值和方差，可以使用不同的方法来计算，例如对数变换、线性平均、动态平均等。
    • 将检测窗口划分为若干个子窗口，每个子窗口的大小和形状可以根据实际应用进行调整。
    • 根据期望的误检概率和背景噪声的统计特性，计算每个子窗口的阈值，例如高斯分布下的阈值可以通过计算高斯分布的概率密度函数得到。
    • 对于每个子窗口，比较信号功率与阈值的大小关系，判断该子窗口内是否存在目标信号。
    • 对于检测到的目标信号，可以进行后续的处理和跟踪，例如目标分离、速度估计等。
    • CFAR算法分为两类：一类是均值类CFAR(CA-CFAR)算法，该类算法应用的前提是假设背景杂波是均匀分布的；另一类是有序统计类CFAR(OS-CFAR)算法，这类算法是为了应对邻域内多目标情况而设计的
  • 我的理解就是一个滑动窗口，其在不断滤波，类似深度学习中的卷积（但不一样，只是都是滤波的感觉）

1.均值类CFAR

image-20231108014322252 $$ \alpha =N(P_{FA}^{-\frac 1N}-1),P_{FA}为虚警概率 $$

1.1 CA-CFAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function [ index, XT ] = cfar_ca( xc, N, pro_N, PAD)

alpha=N.*(PAD.^(-1./N)-1);
index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=(sum(cell_left)+sum(cell_right))./N;
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

1.2 GO-CFAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [ index, XT ] = cfar_go( xc, N, pro_N, PAD)

alpha=N.*(PAD.^(-1./N)-1);
index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=max([mean(cell_left),mean(cell_right)]);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

1.3 SO-CFAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [ index, XT ] = cfar_so( xc, N, pro_N, PAD)

alpha=N.*(PAD.^(-1./N)-1);
index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=min([mean(cell_left),mean(cell_right)]);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

2.统计类CFAR

2.1 OS-CFAR

image-20231108015709762 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [ index, XT ] = cfar_os( xc, N, k, pro_N, PAD)

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

index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    cell_all=cat(2,cell_left,cell_right);
    cell_sort=sort(cell_all);
    
    Z=cell_sort(1,k);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

CFAR算法的MATLAB建模

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
% 函数名：CFAR_1D_Compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 创建人：沈诗艺
  % 日期：2023/11/11
  % 函数功能的简单描述：此函数执行一维CFAR处理，并将处理后得到的阈值与原频谱进行比较，从而得到最终的频谱值
  % 函数特别说明：此函数中对频谱序列边缘均采用补零操作
  % 输入参数：
  %     SF: 经过FFT后的频谱
  %     N_base: 一维CFAR处理的基本单元（两个对称基本单元的和）
  %     N_pro: 一维CFAR处理的保护单元（两个对称保护单元的和）
  %     Pfa：一维CFAR处理的虚警概率
  %     cfar_mode：此函数设置了4种常见的CFAR处理模式：
  %         cfar_mode = 0：执行CA-CFAR
  %         cfar_mode = 1：执行GO-CFAR
  %         cfar_mode = 2：执行SO-CFAR
  %         cfar_mode = 3：执行OS-CFAR
  %     OS_k：若执行OS-CFAR时，需设置取第k个最大值作为阈值，若省略此值，则默认为1
  % 输出参数：
  %     XT: CFAR处理后得到的阈值
  %     SF_cfar：原频谱abs(SF)与阈值XT比较后得到的最终CFAR结果
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 修改人：
  % 日期：
  % 相对上一版程序的修改部分：
  % 修改的内容说明：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [XT, SF_cfar]=CFAR_1D_Compute(SF, N_base, N_pro, Pfa, cfar_mode, OS_k)

    if(~exist('OS_k','var'))
        OS_k = 1;                                                           %如果未出现该变量，则对其进行赋值
    end

    xc_pad=[zeros(1,N_base/2+N_pro/2), abs(SF), zeros(1,N_base/2+N_pro/2)]; %补零操作
    
    alpha=N_base.*(Pfa.^(-1./N_base)-1);
    Z=zeros(1,length(SF));
    XT=zeros(1,length(SF));

    for i=1:1:length(SF)
        cell_left = xc_pad(1,i:i+N_base/2-1);                               %滑动窗左侧的基本单元
        cell_right = xc_pad(1,i+N_base/2+N_pro+1:i+N_base+N_pro);           %滑动窗右侧的基本单元
        if(cfar_mode == 0)                                                  %CA-CFAR
            Z(1,i) = (sum(cell_left)+sum(cell_right))./N_base;
        elseif(cfar_mode == 1)                                              %GO-CFAR
            Z(1,i) = max([mean(cell_left),mean(cell_right)]);
        elseif(cfar_mode == 2)                                              %SO-CFAR
            Z(1,i) = min([mean(cell_left),mean(cell_right)]);
        elseif(cfar_mode == 3)                                              %OS-CFAR
            cell_all = cat(2,cell_left,cell_right);
            cell_sort = sort(cell_all);
            Z(1,i) = cell_sort(1,OS_k);
        end
    end

    XT(1,:)=Z.*alpha;                                                       %其结果是每个点的阈值

    SF_cfar = zeros(1,length(SF));                  
    for i=1:1:length(SF)
        if(abs(SF(i))>XT(i))                                                %频谱>XT，则保留原值
            SF_cfar(1,i) = abs(SF(i));
        else
            SF_cfar(1,i) = 0;                                               %频谱<XT，则对应点频谱置零
        end
    end
end

CFAR算法的Verilog硬件实现

1.基本思路

  • 目前还只是初步优化思路，还没考虑量化以及参数可配置的问题，后续有时间再继续优化

  • 其大致思路是通过时钟计数，并对计数器进行判断，从而对数据进行不同的计算，此外，数据是进到保护单元长度的时候就开始计算

    image-20231221115036640

  • 首先，在时序电路中得想明白一个事，也就是我们在抓取寄存器输出的值时，它总是在时钟上升沿结束之后才更新的，所以我们在时钟上升沿开始去抓取数据时，只能抓取到旧值，只有在再下一个时钟时，才能抓取到上一个时钟更新的值

    image-20231221115107922

  • 同理，可推导出在第12个时钟上升沿到来时，有：

    image-20231221115143827

  • 当cnt>11时，则可以进行加一个新移进来的数以及减一个最旧的数的操作了。对了，同样，还是得寄存左窗口最旧的那个值

    image-20231221115235604

  • 此外，需要设置一个完成标志位，其意味着所有窗口左右两侧的和均算出

    • 简单计算一下：
      • 在cnt == 11时，第一个滑窗左右的和被计算出
      • 此后每个周期都能计算出一此和，这里先不考虑补零问题，那么意味着还需要做100-10-1 = 89次cfar
      • 故当cnt == 11+89(即100时)，所有滑窗求和均算出，此时将finished标志位置高

  • 而最终二者左基本单元和右基本单元的求和计算，会晚一个时钟周期计算出来，道理其实与前面解释的一样，获取时序电路的输出再计算，均会导致晚一个周期，而不是立马被计算出来

2.源代码

  • CA_cfar.v

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    module CA_cfar #(
        parameter  N_reg = 11
    )(
        input clk,
        input rst,
        input start,                //启动信号
        input [6:0] SF,             //输入的FFT后的频谱幅度
        input [4:0] alpha,          //与cfar后Z值相乘的alpha
        output [8:0] XT,            //最终阈值
        output reg finished         //一次cfar完成的标志位
        );

        integer i, j;
        reg [6:0] shift_reg [N_reg - 1 : 0];  //这里的shift_reg右边的位宽实际代表N_reg

        reg [6:0] Ncount;       //时钟计数
        reg [6:0] temp;               //寄存左侧窗口移出的值

        //移位寄存器实现
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                for(i = 0; i < N_reg - 1; i = i + 1) begin
                    shift_reg[i] <= 0;
                end
                temp <= 0;
            end
            else if(start) begin
                for(i = 0; i < N_reg - 1; i = i + 1) begin
                    shift_reg[i] <= 0;
                end  
                temp <= 0; 
            end 
            else begin
                for(j = 0; j < N_reg - 1; j = j + 1) begin
                    shift_reg[j + 1] <= shift_reg[j];
                end
                shift_reg[0] <= SF;
            end
        end

        //时钟计数器
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                Ncount <= 0;
            end
            else if(start) begin
                Ncount <= 0;
            end
            else begin
                Ncount <= Ncount + 1;
            end
        end

        reg valid_flag;            //数据全部移进来的标志信号
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                valid_flag <= 0;
            end
            else begin
                if(Ncount == 11) begin
                    valid_flag <= 1;
                end
                else begin
                    valid_flag <= 0;
                end
            end
        end

        //计算完成标志位
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                finished <= 0;
            end
            else begin
                if(Ncount == 100) begin
                   finished <= 1; 
                end
                else begin
                   finished <= 0;
                end
            end
        end

        //加法求和
        //1.左侧
        reg [8:0] Adder_left;
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                Adder_left <= 9'b0;
            end
            else begin
                if(Ncount == 4) begin
                    Adder_left <= shift_reg[0] + shift_reg[1] + shift_reg[2] + shift_reg[3];
                    temp <= shift_reg[3];
                end
                if(Ncount > 11) begin
                    Adder_left <= Adder_left + shift_reg[7] - temp;
                    temp <= shift_reg[10];
                end
            end
        end   
        //2.右侧
        reg [8:0] Adder_right;
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                Adder_right <= 9'b0;
            end
            else begin 
                if(Ncount == 11) begin
                    Adder_right <= shift_reg[0] + shift_reg[1] + shift_reg[2] + shift_reg[3];
                end
                if(Ncount > 11) begin
                    Adder_right <= Adder_right + shift_reg[0] - shift_reg[4];
                end
            end
        end  

        //左右两侧求和
        reg [9:0] sum_mean;
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                sum_mean <= 10'b0;
            end
            else if(start) begin
                sum_mean <= 10'b0;
            end
            else begin
                sum_mean <= Adder_left + Adder_right;
            end
        end

    endmodule

3.Testbench

  • CA_cfar_tb.v

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    `define CLOCK_PERIOD 10 
    module CA_cfar_tb;

        reg clk;
        reg rst;
        reg start;               
        reg [6:0] SF;            
        reg [4:0] alpha;         
        wire [8:0] XT;         
        wire finished;     

        CA_cfar #(.N_reg(11)) ca_U1(
            .clk(clk),
            .rst(rst),
            .start(start),         
            .SF(SF),            
            .alpha(alpha),         
            .XT(XT),            
            .finished(finished)       
            );

        initial clk = 1;
        always #(`CLOCK_PERIOD/2) clk = ~clk;

        initial begin
            rst = 0;
            start = 0;
            #21 rst = 1;
            #20 start = 1;
            #10 start = 0;
        end

        reg [6:0] dataI[99:0];

        initial begin
            $readmemb("D:\\App_Data_File\\Vivado_data\\Vivado_project\\DSP_design\\CFAR_project\\cfar_test_data.txt",dataI);
        end

        reg [6:0] ncount;
        always @(posedge clk or negedge rst) begin
            if(~rst) begin
                ncount <= 7'b0;
            end
            else if(start) begin
                ncount <= 0;
            end
            else begin
                ncount <= ncount + 1;
            end
        end

        always @(*) begin
            if(~rst) begin
                SF <= 7'b0;
            end
            else begin
                if(ncount < 100) begin
                    SF = dataI[ncount];
                end
                else begin
                    SF = 7'b0;
                end
            end
        end

    endmodule

    • 结果如下：

      image-20231221115618703

      image-20231221115626102

4.更新

  • 目前一维CFAR可以做到CA、GO、SO、OS四种CFAR可选、保护单元和基本单元可调，等项目结束后将有偿提供代码，着急需要代码的可邮箱联系（2024.3.22）

Reference

欢迎来到ssy的世界