利用matlab和FPGA产生FMCW波

了了J

1540人浏览 · 2025-07-01 15:39:47

了了J · 2025-07-01 15:39:47 发布

一、FMCW波

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

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

$S_{m}(t)=Acos[\omega _{c}t+\varphi (t)]$

其中，A是振幅， $[\omega _{c}t+\varphi (t)]$ 是信号的瞬时相位， $\varphi(t)$ 是瞬时相位偏移(相对于 $\omega _{c}t$ )， $\frac{d[[\omega _{c}t+\varphi (t)]]}{dt}$ 是信号的瞬时角频率， $\frac{d\varphi (t)}{dt}$ 是信号的瞬时角频偏(相对于 $\omega _{c}t$ )，除个2 $\pi$ 就从角频率变为频率。

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

$\frac{d\varphi (t)}{dt}\cdot \frac{1}{2\pi }=K_{f}m(t)$

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

$\varphi (t)=2\pi K_{f}\int m(\tau )d\tau$

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

$S_{FM}(t)=Acos[\omega _{c}t+2\pi K_{f}\int m(\tau )d\tau]$

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

$\frac{d\varphi (t)}{dt}\cdot \frac{1}{2\pi }=Kt$

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

$\varphi (t)=2\pi K\int \tau d\tau =2\pi \cdot K\cdot \frac{t^{2}}{2}$

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

$S_(t)=Acos(\omega _{c}t+2\pi\cdot K\cdot \frac{t^{2}}{2} )=Acos[2\pi (f_{c}+\frac{K}{2}t)t]$

将 $2\pi (f_{c}+\frac{K}{2}t)$ 看作是 $\omega$ 的话，可以看出其频率从fc开始以K/2的斜率线性增加。

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

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

如果假设发射信号Si为

$S_{i}=Acos[2\pi (f_{c}+\frac{K}{2}t)t]$

则其频率为

$f_{i}(t)=f_{c}+Kt$

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

$S_{r}=Acos[2\pi (f_{c}+\frac{K}{2}(t-\tau ))(t-\tau )]$

其中 $\tau$ 为时延，则其频率为

$f_{r}(t)=f_{c}+K(t-\tau )$

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

$f_{b}=f_{i}(t)-f_{r}(t)=K\tau$

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

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

$\frac{f_{b}}{\tau }=\frac{B}{T}$

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

$\tau =\frac{2D}{v}$

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

$D=\frac{v\cdot f_{b}\cdot T}{2B}$

也就是说，只要得到频率差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 = 0:dt:Ts-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米\n',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的输出频率有计算公式

$f_{out}=\frac{f_{clk\cdot K}}{2^{N}}$

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

$K=f_{out}\cdot 2^{N}/f_{clk}=25M\cdot 2^{32}/200M=2^{29}=536870912$

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

$K=f_{out}\cdot 2^{N}/f_{clk}=250\cdot 2^{32}/200M\approx 5369$

也就是每个时钟周期频率控制字增加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米\n',d);

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

并且混频后的中频频率差为 $\frac{50}{200000}\cdot 200MHz=50000Hz$ ，与上述matlab产生的数据一样，说明DDS产生的波形无误。

乐奇 Rokid 开放社区

更多推荐

AI智能眼镜功率电感怎么选：顺络SWPA/SPH/WPN/WPG系列与三大工艺平台解析

乐奇 Rokid 开放社区

极刻觅镜 | AI眼镜日报｜Synaptics发布边缘AI处理器、Rokid联手MythLink

乐奇 Rokid 开放社区

Last Call！实时智能 Builder 聚集！RTE Open Day@iRTE 2026，10.23/24，北京

它面向有声剧、视频配音、游戏角色、虚拟主播、智能客服与 Voice Agent，提供自然、可控、实时的语音能力。我们连接开发者、创业者与技术探索者，让新想法在交流中被看见、被验证，也让人与人、人与 AI、人与世界之间的实时互动不断发生。实时人像替换、虚拟换装、风格重渲和虚实融合，可用于直播互动、游戏文娱、智能硬件与电商试穿。还有美团、阿里 ATH、Rokid、斜跃智能、博通集成、Soul、HOLL