调频信号FM的原理与matlab与FPGA实现

平台：matlab r2021b，vivado2023.1

本文知识内容摘自《软件无线电原理和应用》

调频(FM)是载波的瞬时频率随调制信号成线性变化的一种调制方式，音频调频信号的数学表达式可以写为：

$s(t)=A[cos(\omega _{c}t+k_{\Omega }\int_{0}^{t}v_{\Omega (t)dt})]$

Fm频率调制，载波的幅度随着调制波形的幅度变化而变化。

其中 $\omega _{c}$ 为载波频率， $v_{\Omega }t$ 为调制信号， $k_{\Omega }$ 为调制角频率。

下面是FM调制的matlab实现

clc;
clear;
% 设置参数
fs = 312.5e6;       % 采样率
fc = 20e6;          % 载波频率
fm = 1e6;           % 调制信号频率，内调制最大3mhz
t1 = 0:1/fs:2;      % 时间序列,1微秒
t = t1(1:5000);     %RW需要取整数计算出的频率是真实
% 生成调制信号
m = cos(2*pi*fm*t);%正弦波
% m = square(2*pi*fm*t);%方波
% m = sawtooth(2*pi*fm*t, 0.5);%三角波
% m = sawtooth(2*pi*fm*t);% 锯齿波
% 例如，y=x^2;t=1-5;
% Q=cumtrapa(t,y);
% q0=0;
% q1=0.5*(4+1)+0=2.5;
% q2=0.5*(9+3)+2.5=9;
% q3=0.5*(16+9)+9=21.5;
% q4=0.5*(25+16)+21.5=42;
% 计算积分累计积分结果，返回一个向量
integral_term = cumtrapz(t, m);% 生成载波信号
c = cos(2*pi*fc*t);
% FM调制,
kf = 100e6; % 调频系数
k = 2*pi*fc*t;
k1= kf*integral_term;
s = cos(2*pi*fc*t + kf*integral_term);
% 绘制时域波形
figure(1);
subplot(3,1,1);
plot(t*1e6, m);
title('调制信号');
xlabel('时间 (μs)');
ylabel('幅度');subplot(3,1,2);
plot(t*1e6, c);
title('载波信号');
xlabel('时间 (μs)');
ylabel('幅度');subplot(3,1,3);
plot(t*1e6, s);
title('调制后信号');
xlabel('时间 (μs)');
ylabel('幅度');figure(2);
subplot(1,1,1);
plot(t*1e6, c, 'r', 'LineWidth', 2); % 使用红色绘制载波信号，线条宽度为2
title('载波信号');
xlabel('时间 (μs)');
ylabel('幅度');
hold on;plot(t*1e6, s, 'k', 'LineWidth', 2); % 使用黑色绘制调制后信号，线条宽度为2
title('调制后信号');
xlabel('时间 (μs)');
ylabel('幅度');% 绘制频域波形
figure(3);
% 计算频谱
N = length(t);
f = (-fs/2:fs/N:fs/2-fs/N); % 频率向量
M = fftshift(fft(m));
C = fftshift(fft(c));
S = fftshift(fft(s));subplot(3,1,1);
plot(f, abs(M)/N,'b');
title('调制信号频谱');
xlabel('频率 (GHz)');
ylabel('幅度');subplot(3,1,2);
plot(f, abs(C)/N,'g');
title('载波信号频谱');
xlabel('频率 (GHz)');
ylabel('幅度');subplot(3,1,3);
plot(f, abs(S)/N,'r');
title('调制后信号频谱');
xlabel('频率 (GHz)');
ylabel('幅度');

调制波为余弦波时时域和频域图像

当调制信号为方波时。

为锯齿波时

FPGA的实现

当我们的调制波是一个余弦波时。

可以看到我们的调制波形是一个余弦波。载波也是一个余弦波，调制波的频率随着调制波形的积

分变化而变化。其变化规律如下。

余弦波时，积分量在0，pi和2pi时最小，对应着在0时频偏最小，在pi/2时频率与载波相同，在pi时

频偏正向最大。在3*pi/2时又与载波频率相同。在2pi时达到了最小频偏。

在逻辑中有几种产生正余弦波形的方式，基于DDS的波形发生器，基于cordic的波形发生器。这里我们使用cordic来产生我们的载波和调制波。

关于cordic的频率控制字这里说明一下。Cordic是你对其输入一个角度，他给你计算出y（cos，sin）的一个结果。所以我们需要对频率控制字执行一个累加的过程。其中cordic的角度范围表示为（-pi，pi）。

关于输出的频率计算公式为

$f_{o}=\tfrac{phase*fs}{2^{N-2}}$

其中 $f_{o}$ 为输出频率，phase为相位控制字， $f_{s}$ 为采样率。 $2^{N-2}$ 是因为cordic将数据的范围量化到（-pi，pi）。

所以我们需要控制cordic的累加量

$p=p+pi+po$

其中p为频率控制字，pi为载波的频率控制字，po为频偏控制字。

例如我们要载波为fi，最大频偏为fo。假定现在的采样率时钟为fs。根据公式

可以算出载波的频率控制字为

$pi=\frac{{{}2^{14}}*fi}{fs}$

可以算出最大频偏控制字为

$po=\frac{{{}2^{14}}*fo}{fs}$

所以又调制波的幅度最大为16’h4000=16’d16384表示最大为正1v。
所以po与幅度的对应关系为

$k=\frac{{{}2^{14}}*fo}{fs*2^{14}}=\tfrac{fo}{fs}$

所以最大频偏和调制波幅度的关系为

$po=\frac{fo}{fs}*x$

其中x为调制波幅度

逻辑实现现在假定调制波为1mhz，载波为8mhz，最大频偏为2mhz，采样率为512mhz。

插入FPGA代码

`timescale 1ns / 1ps
//
// Company: 
// Engineer: 
// 
// Create Date: 2024/06/06 21:09:44
// Design Name: 
// Module Name: vtf_cordic
// Project Name: 
// Target Devices: 
// Tool Versions: 
// Description: 
// 
// Dependencies: 
// 
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
// 
//module vtf_cordic;reg             aclk;
reg             rst_n;
reg             s_axis_phase_tvalid;
reg  [15 : 0]   s_axis_phase_tdata;
reg  [15 : 0]   s_axis_phase_tdata_0;wire            m_axis_dout_tvalid;
wire [31 : 0]   m_axis_dout_tdata;wire [15:0]     sin ;
wire [15:0]     cos ;
wire [15:0]     sin0 ;
wire [15:0]     cos0 ;cordic_0 u_cordic_0 (.aclk                         (aclk                         ),// input wire aclk.s_axis_phase_tvalid          (s_axis_phase_tvalid          ),// input wire s_axis_phase_tvalid.s_axis_phase_tdata           (s_axis_phase_tdata           ),// input wire [15 : 0] s_axis_phase_tdata.m_axis_dout_tvalid           (m_axis_dout_tvalid           ),// output wire m_axis_dout_tvalid.m_axis_dout_tdata            ({sin,cos}                    )// output wire [31 : 0] m_axis_dout_tdata
);cordic_0 u_cordic_1 (.aclk                         (aclk                         ),// input wire aclk.s_axis_phase_tvalid          (s_axis_phase_tvalid          ),// input wire s_axis_phase_tvalid.s_axis_phase_tdata           (s_axis_phase_tdata_0           ),// input wire [15 : 0] s_axis_phase_tdata.m_axis_dout_tvalid           (                             ),// output wire m_axis_dout_tvalid.m_axis_dout_tdata            ({sin0,cos0}                  )// output wire [31 : 0] m_axis_dout_tdata
);initial
beginaclk    =0;rst_n   =0;#100;rst_n   =1;#100;s_axis_phase_tvalid =1;endreg     [15:0]      wave_add;
reg     [15:0]      phase_tdata;
reg     [15:0]      phase_tdata_0;
//产生一个载波
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginwave_add  <= 16'b0;endelse beginwave_add  <= 16'd32;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginphase_tdata  <= 16'b0;endelse beginphase_tdata  <= phase_tdata + wave_add;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginphase_tdata_0  <= 16'b0;endelse if(phase_tdata >= 16'h0 && phase_tdata <= 16'h4000 )beginphase_tdata_0  <= phase_tdata;endelse if(phase_tdata > 16'h4000 && phase_tdata <= 16'h8000 )beginphase_tdata_0  <= phase_tdata - 16'h4000;endelse if(phase_tdata > 16'h8000 && phase_tdata <= 16'hc000 )beginphase_tdata_0  <= phase_tdata - 16'h8000;endelse if(phase_tdata > 16'hc000 && phase_tdata <= 16'hffff )beginphase_tdata_0  <= phase_tdata - 16'hc000;endelse beginphase_tdata_0  <= phase_tdata;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)begins_axis_phase_tdata  <= 16'b0;endelse begins_axis_phase_tdata  <= 16'he000 + phase_tdata_0;end
end//-------------------------------------------------------------------
reg     [15:0]      wave_add_m;
reg     [15:0]      phase_tdat_m;
reg     [15:0]      phase_tdata_0_m;wire    [15:0]      sinsin={sin[15],sin[15],sin[15],sin[15],sin[15],sin[15],sin[15],sin[15],sin[15:8]};
//产生一个方波
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginwave_add_m  <= 16'b0;endelse beginwave_add_m  <= sinsin + 16'd262;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginphase_tdat_m  <= 16'b0;endelse beginphase_tdat_m  <= phase_tdat_m + wave_add_m;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)beginphase_tdata_0_m  <= 16'b0;endelse if(phase_tdat_m >= 16'h0 && phase_tdat_m <= 16'h4000 )beginphase_tdata_0_m  <= phase_tdat_m;endelse if(phase_tdat_m > 16'h4000 && phase_tdat_m <= 16'h8000 )beginphase_tdata_0_m  <= phase_tdat_m - 16'h4000;endelse if(phase_tdat_m > 16'h8000 && phase_tdat_m <= 16'hc000 )beginphase_tdata_0_m  <= phase_tdat_m - 16'h8000;endelse if(phase_tdat_m > 16'hc000 && phase_tdat_m <= 16'hffff )beginphase_tdata_0_m  <= phase_tdat_m - 16'hc000;endelse beginphase_tdata_0_m  <= phase_tdat_m;end
end
always@(posedge aclk or negedge rst_n)
beginif(rst_n == 1'b0)begins_axis_phase_tdata_0  <= 16'b0;endelse begins_axis_phase_tdata_0  <= 16'he000 + phase_tdata_0_m;end
endalways#0.977 aclk = ~aclk;endmodule

仿真为