既然你能搜到这篇博客,我想你已经大概了解过了脉冲雷达测距的原理,其中最核心的部分就是匹配滤波了,我对于匹配滤波的理解其实就是在计算两个信号的相关性,那么数学上计算两个信号的相关性还是很简单的,下面我就不做过多原理上的介绍了,直接从工程的角度去实现吧。
一、matlab仿真实现
1.1LFM线性调频信号的产生
%%线性调频信号
T=10e-6; %p脉冲持续时间10us
B=30e6; %线性调频信号的频带宽度30MHz
K=B/T; %调频斜率
Fs=2*B;Ts=1/Fs; %采样频率和采样间隔
N=T/Ts; %采样点数
t=linspace(-T/2,T/2,N);
St=exp(1i*pi*K*t.^2); %线性调频信号
subplot(311)
plot(t*1e6,real(St));
xlabel('时间/us');
title('线性调频信号的实部');
subplot(312)
plot(t*1e6,imag(St));
xlabel('时间/us');
title('线性调频信号的虚部');
subplot(313)
freq=linspace(-Fs/2,Fs/2,N);
plot(freq*1e-6,fftshift(abs(fft(St))));
xlabel('频率/MHz');
title('线性调频信号的幅频特性');
figure;
plot(t*1e6,real(St)+imag(St));
xlabel('时间/us');
title('线性调频信号');
上面我们可以看到已经产生了对应的lfm信号,下面我们更具这个信号模拟回波信号的产生。
1.2雷达回波信号的产生
%% 参数初始化
% clc,clear;
C=3e8; %光速
T=50e-6; %脉冲宽度50us
Tr=1e-3; %脉冲重复周期
Fc=30e6; %LFM载频
B=2e6; %频带宽度2MHz
K=B/T; %调频斜率
Fs=40e6;Ts=1/Fs; %采样频率与采样时间
%% 雷达接收窗
Rmin=10e3;Rmax=150e3; %测距范围
R=[100e3]; %目标点的位置,每一个目标相对于雷达的斜距
M=length(R); %目标的个数
RCS=[1]; %雷达截面积,一维数组
Rwid=Rmax-Rmin; %最大测距长度
Twid=2*Rwid/C; %回波窗的长度
Nwid=ceil(Twid/Ts); %采样窗内的采样点数
t=linspace(2*Rmin/C,2*Rmax/C,Nwid); %回波窗
%% 产生回波
td=ones(M,1)*t-2*R'/C*ones(1,Nwid); %回波时间序列
Srt=real(RCS*(exp(2j*pi*(Fc*td+1/2*K*td.^2)).*(abs(td)<T/2)));%从点目标来的回波
figure;
plot(Srt);
%axis([0 35000 -1 1]);
xlabel('时间/us');
title('回波信号');
1.3雷达回波信号的量化
因为要在fpga种去实现,所以要将该回波信号量化成.coe文件
%% 参数初始化
clc,clear;
C=3e8; %光速
T=50e-6; %脉冲宽度50us
Tr=1e-3; %脉冲重复周期
Fc=30e6; %LFM载频
B=2e6; %频带宽度2MHz
K=B/T; %调频斜率
Fs=40e6;Ts=1/Fs; %采样频率与采样时间
%% 雷达接收窗
Rmin=10e3;Rmax=150e3; %测距范围
R=[100e3]; %目标点的位置,每一个目标相对于雷达的斜距
M=length(R); %目标的个数
RCS=[1]; %雷达截面积,一维数组
Rwid=Rmax-Rmin; %最大测距长度
Twid=2*Rwid/C; %回波窗的长度
Nwid=ceil(Twid/Ts); %采样窗内的采样点数
t=linspace(2*Rmin/C,2*Rmax/C,Nwid); %回波窗
%% 产生回波
td=ones(M,1)*t-2*R'/C*ones(1,Nwid); %回波时间序列
Srt=RCS*(exp(2j*pi*(Fc*td+1/2*K*td.^2)).*(abs(td)<T/2));%从点目标来的回波
Srt_awgn = awgn(real(Srt),60); %回波加60db噪声
Wgn = Srt_awgn - Srt; %加噪回波-回波=噪声
%% 生成回波接收数据.coe文件
width_sig=16; %信号的数据宽度(量化位数)
width_noi= 6; %噪声的数据宽度(量化位数)
Srt_quan= floor(real(Srt) * (2^(width_sig-1)-1))+floor(real(Wgn) * (2^(width_noi-1)-1));%量化
fid = fopen('C:\Users\Lenovo\Desktop\LFM\coe\data.txt','w');%生成.txt文件
fprintf(fid,'MEMORY_INITIALIZATION_RADIX=10;\n');
fprintf(fid,'MEMORY_INITIALIZATION_VECTOR=\n');
fprintf(fid,'%d,\n',Srt_quan);%将采样数据输入文件中
fclose(fid);%关闭文件
在这里我先量化成了.txt文件,主要是先查看量化好的文件是否满足Xilinx所要求的文件格式,然后直接把文件的名字重命名成.coe就好了。
1.4DDC数字下变频
在这里我们使用matlab的fdatool工具产生滤波器系数,并将该系数放入到Num数组中。
%% 数字下变频
%频移
n = 1:length(Srt);
cos_unit = [1,0,-1,0];
cos_osc =cos_unit(mod(n,4)+1);
sin_unit = [0,1,0,-1];
sin_osc =sin_unit(mod(n,4)+1);
DDC_out_i1 = Srt.*cos_osc;
DDC_out_q1 = Srt.*sin_osc;
figure(1);
subplot(211);
plot(t*1e6,DDC_out_i1);
axis([600 700 -2 2]);
subplot(212);
plot(t*1e6,DDC_out_q1);
axis([600 700 -2 2]);
%滤波
DDC_out_i = conv(Num,DDC_out_i1);
DDC_out_q = conv(Num,DDC_out_q1);
figure(2);
subplot(211);
plot(DDC_out_i);
axis([22e3 26e3 -1 1]);
subplot(212);
plot(DDC_out_q);
axis([22e3 26e3 -1 1]);
%抽取
deciN = 10;
DDC_out_i_deci = DDC_out_i(1:deciN:length(DDC_out_i));
DDC_out_q_deci = DDC_out_q(1:deciN:length(DDC_out_q));
figure(3);
subplot(211);
plot(DDC_out_i_deci);
axis([22e2 26e2 -1 1]);
subplot(212);
plot(DDC_out_q_deci);
axis([22e2 26e2 -1 1]);
我们可以得到DDC数字下变频信号为:
对下变频后的信号进行滤波,相当于取了一个包络:
对滤波后的信号进行抽取,目的是为了减少数据量,在这里是10倍抽取,比如原来有100个数据点,抽取完以后就剩10个数据点,得到减少了信号处理的数量。
1.5匹配滤波
%% 脉冲压缩
t_base = -T/2:1/Fs:T/2; %时间序列
h1 = exp(1j*K*pi*t_base.^2);%基带信号
h = conj(fliplr(h1)); %基带信号的共轭反转
hi = real(h);%匹配系数实部
hq = imag(h);%匹配系数虚部
width = 16;%数据宽度(量化位数)
hiw = floor(hi' * (2^(width-1)-1) .*hamming(length(hi)));%系数实部加窗后量化
hqw = floor(hq' * (2^(width-1)-1) .*hamming(length(hq)));%系数虚部加窗后量化
hiw_deci = hiw(1:deciN:length(hiw));
hqw_deci = hqw(1:deciN:length(hqw));
PC_hisi = conv(hiw_deci,DDC_out_i_deci);
PC_hisq = conv(hiw_deci,DDC_out_q_deci);
PC_hqsi = conv(hqw_deci,DDC_out_i_deci);
PC_hqsq = conv(hqw_deci,DDC_out_q_deci);
figure(4);
subplot(411);
plot(PC_hisi);
axis([2200 2800 -200000 200000]);
subplot(412);
plot(PC_hisq);
axis([2200 2800 -200000 200000]);
subplot(413);
plot(PC_hqsi);
axis([2200 2800 -200000 200000]);
subplot(414);
plot(PC_hqsq);
axis([2200 2800 -200000 200000]);
PC_i = PC_hisi - PC_hqsq;
PC_q = PC_hisq + PC_hqsi;
PC_final = sqrt(PC_i.^2 + PC_q.^2) ;
% PC_result = PC_i + 1j*PC_q;
% PC_final = db(abs(PC_result));
dis=linspace(Rmin,Rmax,length(PC_final)); %回波窗
figure(5);
plot(dis,PC_final);
axis([Rmin Rmax 0 200000]);
加入汉明窗做旁瓣压制,且对匹配滤波器系数进行量化后的输出:
目标检测结果:
我们可以看到目标大概现在10*10e4的位置处,就是100KM,与回波的设置相同。
1.6匹配滤波系数量化
%% 参数初始化
clc,clear;
C=3e8; %光速
T=50e-6; %脉冲宽度50us
Tr=1e-3; %脉冲重复周期
Fc=30e6; %LFM载频
B=2e6; %频带宽度2MHz
K=B/T; %调频斜率
Fs=4e6;Ts=1/Fs; %采样频率与采样时间
%% 生成本地匹配滤波器
deciN = 10;
width = 16;%数据宽度(量化位数)
t = -T/2:1/Fs:T/2; %时间序列
h = exp(-1j*K*pi*t.^2);%信号基带复共轭
hi = real(h);%匹配系数实部
hq = imag(h);%匹配系数虚部
hiw = floor(hi' * (2^(width-1)-1).*hamming(length(hi)));%系数实部加窗后量化
hqw = floor(hq' * (2^(width-1)-1).*hamming(length(hq)));%系数虚部加窗后量化
hiw_deci = hiw(1:deciN:length(hiw));
hqw_deci = hiw(1:deciN:length(hqw));
%% 生成匹配系数实部.coe文件
fid_i = fopen('C:\Users\Lenovo\Desktop\LFM\coe\GeneratePc_I.txt','w');%生成.txt文件
fprintf(fid_i,'Radix = 10; \n');
fprintf(fid_i,'Coefficient_Width = 16; \n');
fprintf(fid_i,'CoefData = \n');
fprintf(fid_i,'%d,\n',hiw_deci);%将采样数据输入文件中
fclose all;%关闭文件
%% 生成匹配系数虚部.coe文件
fid_q = fopen('C:\Users\Lenovo\Desktop\LFM\coe\GeneratePc_Q.txt','w');%生成.txt文件
fprintf(fid_q,'Radix = 10; \n');
fprintf(fid_q,'Coefficient_Width = 16; \n');
fprintf(fid_q,'CoefData = \n');
fprintf(fid_q,'%d,\n',hqw_deci);%将采样数据输入文件中
fclose all;%关闭文件
和雷达回波信号的量化一样,都是为了产生让Xilinx所要求的数据格式,注意量化完成后要将文佳后缀名改成.coe。
二、FPGA实现
2.1顶层文件
`timescale 1ns / 1ps
module top_test
(
input wire clk_100M,
input wire rst_n,
output wire [47:0]PC_final //[54:0]
);
wire [15:0] data; //ram中输出数据
wire [15:0] DDC_out_i1; //i路未滤波信号
wire [15:0] DDC_out_q1; //q路未滤波信号
wire [34:0]DDC_out_i; //i路滤波信号
wire [34:0]DDC_out_q; //q路滤波信号
wire valid_i; //i路抽取信号
wire valid_q; //q路抽取信号
wire [53:0]PC_hisi;
wire [53:0]PC_hisq;
wire [53:0]PC_hqsq;
wire [53:0]PC_hqsi;
data_gen data_gen_inst
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.data_out(data)
);
mix_fre mix_fre_inst
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.data_in(data),
.DDC_out_i1(DDC_out_i1),
.DDC_out_q1(DDC_out_q1)
);
ddc_fir ddc_fir_inst_i
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.dina(DDC_out_i1),
.dout(DDC_out_i)
);//I路低通滤波
ddc_fir ddc_fir_inst_q
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.dina(DDC_out_q1),
.dout(DDC_out_q)
);//Q路低通滤波
chouqu chouqu_inst_i
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.valid(valid_i)
);
chouqu chouqu_inst_q
(
.clk_100M(clk_100M),
.rst_n(rst_n),
.valid(valid_q)
);
PulseCompress_hi PulseCompress_hisi
(
.clk_100M(clk_100M), //系统时钟100MHz
.s_axis_data_tvalid(valid_i),//数据有效电平
.din(DDC_out_i), //数据输入
.dout(PC_hisi) //数据输出
);
PulseCompress_hi PulseCompress_hisq
(
.clk_100M(clk_100M), //系统时钟100MHz
.s_axis_data_tvalid(valid_i),//数据有效电平
.din(DDC_out_i), //数据输入
.dout(PC_hisq) //数据输出
);
PulseCompress_hq PulseCompress_hqsi
(
.clk_100M(clk_100M), //系统时钟100MHz
.s_axis_data_tvalid(valid_q),//数据有效电平
.din(DDC_out_q), //数据输入
.dout(PC_hqsi) //数据输出
);
PulseCompress_hq PulseCompress_hqsq
(
.clk_100M(clk_100M), //系统时钟100MHz
.s_axis_data_tvalid(valid_q),//数据有效电平
.din(DDC_out_q), //数据输入
.dout(PC_hqsq) //数据输出
);
calue calue_inst
(
.PC_hisi(PC_hisi),
.PC_hisq(PC_hisq),
.PC_hqsq(PC_hqsq),
.PC_hqsi(PC_hqsi),
.PC_final(PC_final)
);
endmodule
在顶层中我们可以看到,这里包含了一个data_gen模块,这个模块是将量化好的雷达回波信号保存在这里,然后再读出;还包含一个mix_fre模块,这个模块是将雷达回波信号进行DDC下变频处理,输出I路和Q路的信号;还包含两个ddc_fir模块,这个模块是将I路和Q路信号进行滤波处理,输出DDC_out_i和DDC_out_q信号;还包含两个chouqu模块,这个模块是将滤波后的I路和路信号进行抽取,输出valid_i和valid_q;还包含两个PulseCompress_hi模块,这个是将I路信号进行匹配滤波;还包含两个PulseCompress_hq模块,这个是将Q路信号进行匹配滤波处理;最好包含一个calue模块,这个模块是按照复数的乘法规则进行各路信号的相加减运算得到脉压后的实部和虚部信号最后得到目标信号。
2.2data_gen模块
`timescale 1ns / 1ps
module data_gen
(
input wire clk_100M,
input wire rst_n,
output wire [15:0]data_out
);
//wire [15:0] dina; //输入信号
reg [15:0] addr; //地址初始化
//wire signed [15:0] dout;//读出回波数据
//按时钟读出回波信号数据
always @ (posedge clk_100M or negedge rst_n)
begin
if(rst_n == 1'b0)
addr <= 16'd0;
else if(addr == 16'd37334) //注意此处地址上限值是数据深度
addr <= 16'd0;
else
addr <= addr + 16'b1;
end
BRAM BRAM_inst (
.clka(clk_100M), // input wire clka
.ena(1'b1), // input wire ena
.wea(1'b0), // input wire [0 : 0] wea
.addra(addr), // input wire [15 : 0] addra
.dina(), // input wire [15 : 0] dina
.douta(data_out) // output wire [15 : 0] douta
);
endmodule
这个里面例化了一个ram,里面保存到数据就是雷达回波信号数据。
上图就是FPGA所产生的雷达回波信号。
2.3mix_fre模块
`timescale 1ns / 1ps
module mix_fre
(
input wire clk_100M,
input wire rst_n,
input wire [15:0]data_in,
output reg [15:0]DDC_out_i1,
output reg [15:0]DDC_out_q1
);
/*混频:由于采用的是最优采样频率(4f_0/(2m+1)),又根据取的m值为1,
因此I、Q路的本振信号分别是[1,0,-1,0]、[0,1,0,-1]的周期序列,
若不满足最优采样频率可添加DDS进行混频*/
reg [1:0] count ; //计数器用以计数状态
//reg signed [15:0] DDC_out_i1=0; //混频后的I路信号
//reg signed [15:0] DDC_out_q1=0; //混频后的Q路信号
wire signed [15:0] AD_in_p; //回波原信号
wire signed [15:0] AD_in_n; //回波原信号相反数
assign AD_in_p = data_in;
assign AD_in_n = ~data_in[15:0] + 16'd1;//!硬件语言中的取相反数操作
always@(posedge clk_100M or negedge rst_n)
begin
if(rst_n == 1'b0)
count <= 2'b0;
else if(count == 2'd3)
count <= 2'b0;
else
count <= count + 1'b1;
case(count)
2'd0:begin
DDC_out_i1 <= AD_in_p;
DDC_out_q1 <= 0;
end
2'd1:begin
DDC_out_i1 <= 0;
DDC_out_q1 <= AD_in_p;
end
2'd2:begin
DDC_out_i1 <= AD_in_n;
DDC_out_q1 <= 0;
end
2'd3:begin
DDC_out_i1 <= 0;
DDC_out_q1 <= AD_in_n;
end
endcase
end
endmodule
2.4ddc_fir模块
`timescale 1ns / 1ps
module ddc_fir
(
input wire clk_100M,//系统时钟100MHz
input wire rst_n ,//复位按钮(低电平有效)
input wire signed [15:0] dina ,//回波数据信号输入
output wire signed [34:0] dout //滤波输出
);
//滤波器初始化
reg s_axis_data_tvalid;
//wire s_axis_data_tready;
//wire m_axis_data_tvalid;
wire signed [39:0] m_axis_data_tdata;
//复位时拉高数据有效电平
always @ (posedge clk_100M or negedge rst_n)
begin
if(rst_n == 1'b0)
s_axis_data_tvalid <= 1'b0;
else if(rst_n == 1'b1)
s_axis_data_tvalid <= 1'b1;
else
s_axis_data_tvalid <= s_axis_data_tvalid;
end
//FIR的IP核例化
fir fir_inst
(
.aclk(clk_100M), // input wire aclk
.s_axis_data_tvalid(s_axis_data_tvalid), // input wire s_axis_data_tvalid
.s_axis_data_tready(), // output wire s_axis_data_tready
.s_axis_data_tdata(dina), // input wire [15 : 0] s_axis_data_tdata
.m_axis_data_tvalid(), // output wire m_axis_data_tvalid
.m_axis_data_tdata(m_axis_data_tdata) // output wire [39 : 0] m_axis_data_tdata
);
assign dout = m_axis_data_tdata[34:0];//截取输出信号的低35位
endmodule
2.5chouqu模块
`timescale 1ns / 1ps
module chouqu
(
input wire clk_100M,
input wire rst_n,
output reg valid
);
reg [3:0]cnt;
//10分频,相当于10倍抽取
always @ (posedge clk_100M or negedge rst_n)
begin
if(rst_n == 1'b0)
cnt <= 4'b0;
else if((rst_n == 1'b1) && (cnt == 4'd9))
cnt <= 4'b0;
else
cnt <= cnt + 4'b1;
end
always @ (posedge clk_100M or negedge rst_n)
begin
if(rst_n == 1'b0)
valid <= 1'b0;
else if((cnt == 4'd9) && (rst_n == 1'b1))
valid <= 1'b1;
else
valid <= 1'b0;
end
endmodule
在这个模块中没10个时钟单位产生一次valid信号,这个valid用于控制后面匹配滤波器的使能信号。(在这FPGA产生的信号就不展示了)
2.6PulseCompress_hi和PulseCompress_hq模块
这两个模块是匹配滤波模块
`timescale 1ns / 1ps
module PulseCompress_hi(
input clk_100M, //系统时钟100MHz
input s_axis_data_tvalid,//数据有效电平
input wire signed [34:0] din, //数据输入
output wire signed [53:0] dout //数据输出
);
//FIR引脚初始化
//wire s_axis_data_tready;
//wire m_axis_data_tvalid;
wire signed [55:0] m_axis_data_tdata ;
//FIR例化
compress_i compress_i_inst
(
.aclk(clk_100M), // input wire aclk
.s_axis_data_tvalid(s_axis_data_tvalid), // input wire s_axis_data_tvalid
.s_axis_data_tready(), // output wire s_axis_data_tready
.s_axis_data_tdata({{5{din[34]}},din}), // input wire [39 : 0] s_axis_data_tdata
.m_axis_data_tvalid(), // output wire m_axis_data_tvalid
.m_axis_data_tdata(m_axis_data_tdata) // output wire [55 : 0] m_axis_data_tdata
);
assign dout = m_axis_data_tdata[53:0];//截取FIR输出的低54位
endmodule
`timescale 1ns / 1ps
module PulseCompress_hq(
input clk_100M, //系统时钟100MHz
input s_axis_data_tvalid,//数据有效电平
input wire signed [34:0] din, //数据输入
output wire signed [53:0] dout //数据输出
);
//FIR引脚初始化
//wire s_axis_data_tready;
//wire m_axis_data_tvalid;
wire signed [55:0] m_axis_data_tdata ;
//FIR例化
compress_q compres_q_inst
(
.aclk(clk_100M), // input wire aclk
.s_axis_data_tvalid(s_axis_data_tvalid), // input wire s_axis_data_tvalid
.s_axis_data_tready(), // output wire s_axis_data_tready
.s_axis_data_tdata({{5{din[34]}},din}), // input wire [39 : 0] s_axis_data_tdata
.m_axis_data_tvalid(), // output wire m_axis_data_tvalid
.m_axis_data_tdata(m_axis_data_tdata) // output wire [55 : 0] m_axis_data_tdata
);
assign dout = m_axis_data_tdata[53:0];//截取FIR输出的低54位
endmodule
匹配滤波后的输出如上图所示,但是感觉不太理想,我想应该是量化的时候处理问题。
2.6calue模块
`timescale 1ns / 1ps
module calue
(
input wire [53:0]PC_hisi,
input wire [53:0]PC_hisq,
input wire [53:0]PC_hqsq,
input wire [53:0]PC_hqsi,
output wire [47:0]PC_final //[54:0]
);
//两路脉冲压缩后的信号
wire signed [54:0] PC_i;//实部信号
wire signed [54:0] PC_q;//虚部信号
//wire [54:0] PC_i_reg;//实部信号绝对值
//wire [54:0] PC_q_reg;//虚部信号绝对值
//按照复数的乘法规则进行各路信号的相加减运算得到脉压后的实部和虚部信号
assign PC_i = {PC_hisi[53],PC_hisi} - {PC_hqsq[53],PC_hqsq};//实部信号
assign PC_q = {PC_hisq[53],PC_hisq} + {PC_hqsi[53],PC_hqsi};//虚部信号
assign PC_final = $signed(PC_i[54:31]) * $signed(PC_i[54:31]) + $signed(PC_q[54:31]) * $signed(PC_q[54:31]);
//assign PC_i_reg = PC_i[54] ? ~PC_i + 1 : PC_i;
//assign PC_q_reg = PC_q[54] ? ~PC_q + 1 : PC_q;
//assign PC_final = PC_i_reg + PC_q_reg;
endmodule
所以计算出的近距离为:
三、误差分析及结语
3.1误差分析
从FPGA的仿真中我们可以看到,最后的输出是有一定误差的,这个误差从哪里来呢?我认为第一点是在最后的计算模块出现的,在这里采用用来一种近似的算法,实际的算法是平方然后开根号,但是这样有会使工程更加复杂,因为FPGA进行开根号运算其实是特别消耗逻辑资源的;第二点,可以看到目标回波的左右多了两个峰,这个波形感觉就像目标被假目标干扰了,我们可以在后面添加一个门限检测模块,以便于提高检测的准确性;第三点,可能是匹配滤波的系数量化出现了问题,我们可以对比matlab的仿真图和fpga的实际波形,可以看出区别还是很大的,加什么类型的窗(这里是汉明窗)还有量化的位数,或多或少都会对结果产生不同的影响。下面给出不同的处理方式matlab的不同结果:
上图是没有加窗的结果。
上图是更换了量化系数后的结果。
但然还有很多不同的组合,这里只是说明可能出现误差的步骤。同时我们知道FPGA无法处理浮点数只能处理定点数所以要怎么量化,量化多少值得更进一步的研究。
3.2结束
首先感谢大家可以看到这里,同时由于本人的能力有限,希望大佬指正。