利用matlab和FPGA产生FMCW波

一、FMCW波

调频连续波(frequency modulated continuous wave, FMCW)，常用于雷达中进行距离检测，其所用到的信号是一种频率随调制信号线性增长的调频波，还可用于电缆的故障定位。本文基于电缆的故障定位这个方向，简单介绍FMCW波的产生，也是对最近学习的一个记录和总结。

由角调信号出发，其一般表达式为

其中，A是振幅， 是信号的瞬时相位， 是瞬时相位偏移(相对于  )， 是信号的瞬时角频率， 是信号的瞬时角频偏(相对于  )，除个2就从角频率变为频率。

对于调频信号，就是瞬时频率偏移随调制信号m(t)成比例变化，即

其中Kf为调频灵敏度，单位是rad/s*V，此时瞬时相位偏移就是

代入角调信号的一般表达式可得调频信号的一般表达式为

上面提到FMCW用的信号是频率随调制信号线性增长的调频波，既然是调频波，频率又是线性增长，则可以理解为瞬时频率偏移线性增长，即

其中K是斜率，则此时瞬时相位偏移就是

于是便可得到线性调频连续波的表达式为

将看作是  的话，可以看出其频率从fc开始以K/2的斜率线性增加。

对于FMCW雷达测距更加详细的原理请看：

调频连续波（FMCW）雷达测距、测速详解-CSDN博客

如果假设发射信号Si为

则其频率为

发射信号进入电缆中碰到故障位置会发生反射，则接收信号Sr为 

其中为时延，则其频率为

若将接收信号和发射信号进行混频，滤波，可以得到一个它们之间的一个频率差，为

也就是说发射信号和接收信号混合后会产生一个频率为的信号。也可以将Si与Sr的三角表达式进行相乘(混频)，也能得到一个频率为的信号。

又由于对于线性调频信号，

其中fb是频率差，为时延，B为带宽，T为调制周期，并且

其中D是故障点距发射点的距离，由于反射，故距离为2D，v是信号在电缆中的传播速度，则结合上述式子，可以得到故障点的距离为

也就是说，只要得到频率差fb和传播速度v就能得到故障点的距离D。

二、matlab产生FMCW波

设置的参数有调制周期Ts为1ms，带宽Bw为50MHz，起始频率fc为25MHz，采样频率fs为200MHz，则斜率tao为Bw/Ts。

clear;clc;

Ts = 0.001;    %调制周期1ms

Bw = 50e6;     %带宽50MHZ

fc=25e6;       %起始频率25MHz

tao = Bw/Ts;   % 斜率 50HZ/ns     250HZ/5ns    fc+tao*t =>(积分) fc*t+(tao/2)*t^2  

fs = 4*Bw;     %采样频率

dt = 1/fs;     %采样间隔             

N = Ts/dt;     %一个周期的采样点数       

t = 0Ts-dt;  %时间序列

Am=1;                       % 信号幅度 

SFM=Am*cos(2*pi*(fc+(tao/2)*t).*t); %发射信号

figure(1);                  % 绘制第一幅图

subplot(211);               % 窗口分割，将一幅图分割成2*2的

plot(t,SFM,'Linewidth',2);   % 时间t为横坐标，基带信号SFM为纵坐标绘图，线宽为2

xlabel('t/时间');           % 横坐标标注

ylabel('幅度');             % 纵坐标标注

title('发射信号');          % 图标题标注

axis([10e-4*Ts*100,10e-4*Ts*101,-1.1,1.1]);     % 横纵坐标范围设置

line([10e-4*Ts*100,10e-4*Ts*101],[0,0],'color','b','Linewidth',2);% 绘制一条从(0,0)到(0.1,0)的蓝色实线，线宽为2

subplot(212);

[mf,msf]=T2F(t,SFM);              % 傅里叶变换，得到纵坐标频谱和横坐标频率

subplot(212);

plot(mf,abs(msf),'Linewidth',1);  % 画出SFM信号频谱    频谱以fm为间隔

title('发射信号的频谱');            % 图标题标注

xlabel('f/Hz');                     % 横坐标标注

ylabel('幅度/H(f)');                % 纵坐标标注

t1=1e-6; %时延 1us     d=c/2t1=150m

SFM1=Am*cos(2*pi*(fc+(tao/2)*(t-t1)).*(t-t1)); %接收信号

Fb=SFM.*SFM1;  %混频

[b,a]=butter(6,fc/(fs/2));        %巴特沃斯滤波器参数

Xb=filter(b,a,Fb);                %滤波      tao*t1的中频频率

figure(2);                          % 绘制第二幅图

subplot(211);                       % 窗口分割，将一幅图分割成2*1的

plot(t,Xb,'Linewidth',2);         

title('接收混频信号');                % 图标题标注

xlabel('t/时间');                   % 横坐标标注

ylabel('幅度');                     % 纵坐标标注

axis([0,Ts*0.5,-1.1,1.1]);     % 横纵坐标范围设置

[mf2,msf2]=T2F(t,Xb);              % 傅里叶变换，得到纵坐标频谱和横坐标频率

subplot(212);

plot(mf2,real(msf2),'Linewidth',1);  % 画出AM信号频谱    频谱以fm为间隔

title('混频信号的频谱');            % 图标题标注

xlabel('f/Hz');                     % 横坐标标注

ylabel('幅度/H(f)');                % 纵坐标标注

d=(3e8*Ts*tao*t1)/(2*Bw);            %距离计算

fprintf('距离为：%d米 ',d);

function [f,sf]= T2F(t,st)      % FFT

% dt = t(2)-t(1);

T=t(end);                       % 输入信号的时间最大值为T

df = 1/T;                       % dt=1/fs; 时间采样间隔，采样频率的倒数;

% N=T/dt;  采样点个数，总时长除以采样间隔

% 两式联合推导 df = 1/T 

N = length(st);                 % 输入信号时间的长度为采样点数

f=-N/2*df : df : N/2 * df-df;   % 频率分布

sf = fft(st);                   % 做FFT

sf = T/N * fftshift(sf);        % 最后输出，将0-fs频谱搬移到-fs/2-fs/2频谱

发射信号的波形图如下，只看了一小部分，实际上观看同样长度，不同时间频率是变化的。

发射信号频谱图如下，可以看到其频率范围为25MHz至75MHz

混频滤波之后的信号如下，从其频谱图可以看到产生了一个50000Hz的信号，由于时延t1设置为1us，频率差tao*t1也是为50000

假设信号在电缆中以光速传播，即v=300000000m/s，则可以计算得到距离D为150m。并且假设雷达信号在空气中以光速传播，则1us后接收到信号，实际是雷达信号传播了300m，即在150m处反射。可以说明上述的公式推导正确。

三、FPGA产生FMCW波

要在FPGA中去产生FMCW波需要用到DDS IP核。利用其频率控制字去控制DDS产生频率线性增长的波。本文利用vivado软件，Verilog语言进行仿真产生波形数据，再将数据导入matlab中进行校验。选取的板块时钟频率为200MHz，即5ns一个数据。

下面是DDS IP核的简单配置，更加详细的配置可以参考其他文章。第一个界面中将参数选择设置为硬件参数，在相位位宽处可以设置需要的位宽，这里是32位。位宽数用于计数频率控制字。

第二个页面将相位的增加形式设置为Streaming，可以在代码中去控制输入的相位数据改变频率，适合产生调频信号。

后续都保持默认就行。

对于DDS的输出频率有计算公式

其中fclk为时钟频率，k为相位控制字，也就是输入的相位数据，N为相位位宽。按照上述matlab中的参数，要想fc为25MHz，则

也就是频率控制字起始为536870912。并且由于调制周期是1ms，时钟频率为200MHz，即1ms有200000个数据，带宽为50MHz，即每个数据增加250Hz，也就是步长为250Hz，则

也就是每个时钟周期频率控制字增加5369。

//产生连续的线性调频信号，200MHz时钟，1ms调制周期，25MHz起始频率，带宽50MHz，步长(一个时钟周期5ns增加)250Hz

//32位，fout=fclk*K/2^32,K位相位控制字，则起始频率K=536870912，步长K=5369

module Gen_LFM

#( 

parameter CNT_TAO_MAX = 32'd199999      //1ms计数值

)

(

input    wire  sys_clk,

input    wire  sys_rst_n,

input    wire  ctl_tao,                //开始信号

output   wire  [16-1:0]data_out

);

reg     [18-1:0]cnt_tao;                  //1ms计数器

reg    tao_en;                            

reg    tao_en_d1;

reg     [32-1:0]dds_pinc_in;

wire    [32-1:0]s_axis_phase_tdata;

parameter   FREQ_START = 32'd536870912;  //起始频率为25MHz

parameter   FREQ_STEP = 32'd5369;       //步进频率为250Hz

wire   [32-1:0] data;

wire  [16-1:0]sin_imag;

wire  [16-1:0]cos_real;

always @(posedge sys_clk or negedge sys_rst_n) begin

if (sys_rst_n == 1'b0) begin

tao_en <= 1'b0;

end

else begin

if (ctl_tao == 1'b1) begin

tao_en <= 1'b1;

end

end

end

always @(posedge sys_clk or negedge sys_rst_n) begin

if (sys_rst_n == 1'b0) begin

tao_en_d1 <= 1'b0;

end

else begin

tao_en_d1 <= tao_en;  // 延迟一个时钟周期

end

end

always @(posedge sys_clk or negedge sys_rst_n) begin

if (sys_rst_n == 1'b0) begin

cnt_tao <= 18'd0;

end

else if (tao_en == 1'b1) begin

if (cnt_tao == CNT_TAO_MAX) begin

cnt_tao <= 18'd0;

tao_en <= 1'b0;

end

else begin

cnt_tao <= cnt_tao + 1'b1;

end

end

end

//频率控制字

always @(posedge sys_clk or negedge sys_rst_n) begin

if (sys_rst_n == 1'b0) begin

dds_pinc_in <= FREQ_START;

end

else begin

if (tao_en == 1'b1) begin

dds_pinc_in <= dds_pinc_in + FREQ_STEP;

end 

else if (tao_en == 1'b0 && tao_en_d1 == 1'b1) begin

dds_pinc_in <= FREQ_START;

end

end 

end

assign s_axis_phase_tdata=dds_pinc_in; 

dds dds_inst (

.aclk(sys_clk),                                // input wire aclk

.aclken(tao_en),                            // input wire aclken

.s_axis_phase_tvalid(tao_en),  // input wire s_axis_phase_tvalid

.s_axis_phase_tdata(s_axis_phase_tdata),    // input wire [63 : 0] s_axis_phase_tdata

.m_axis_data_tvalid(),    // output wire m_axis_data_tvalid

.m_axis_data_tdata(data),      // output wire [15 : 0] m_axis_data_tdata

.m_axis_phase_tvalid(),  // output wire m_axis_phase_tvalid

.m_axis_phase_tdata()    // output wire [31 : 0] m_axis_phase_tdata

);

assign sin_imag = data[31:16];

assign cos_real = data[15:0]+4096 ;

//将数据导出为txt

integer real_out_file;

initial

begin

real_out_file = $fopen("D:/FMCW/FM/cos_real.txt","w");//定义需要保存的路径和文件名

end

always @(posedge sys_clk or  negedge sys_rst_n) begin

if(tao_en)begin

$fdisplay(real_out_file,"%b",cos_real);//将数据保存为二进制

end

else if (tao_en == 1'b0 && tao_en_d1 == 1'b1) begin

$fclose(real_out_file);

end

end

endmodule

在vivado的仿真图中可以看到输出的频率是逐渐增加的，并且在1ms后频率又回到25MHz。

将导出的数据用matlab读取

clear;clc;

Bw = 50e6;

Ts = 0.001; 

c=3e8;

fileID = fopen('cos_real.txt', 'r');

% 读取整个文件内容为字符串

file_content = fscanf(fileID, '%c');

% 关闭文件

fclose(fileID);

% 移除所有非二进制字符（如换行符、空格等）

file_content = regexprep(file_content, '[^01]', '');

% 设置二进制位数（例如8位表示一个字节）

bits_per_number = 16;

% 计算可转换的完整二进制数数量

num_numbers = floor(length(file_content) / bits_per_number);

% 初始化结果数组

decimal_values = zeros(num_numbers, 1);

% 逐块转换

for i = 1:num_numbers

% 提取固定位数的二进制字符串

binary_str = file_content((i-1)*bits_per_number+1 : i*bits_per_number);

% 转换为十进制

decimal_values(i) = bin2dec(binary_str);

end

SFM=(decimal_values-4096)./2048;

SFM1 = [SFM(end-199:end); SFM(1:end-200)];   % 延迟1us(5ns一个数据，即200)

figure(1)

plot(abs(fft(SFM)));

title('发射信号的频谱');            

figure(2)

Fb=SFM.*SFM1;

plot(real(fft(Fb)));

title('混频信号的频谱');           

d=(50000*Ts*c)/(2*Bw);

fprintf('距离为：%d米 ',d);

 

 

采样频率为时钟频率200MHz，从频谱图可以看出产生的信号频率在25MHz至75MHz。

并且混频后的中频频率差为，与上述matlab产生的数据一样，说明DDS产生的波形无误。