基于全相位FFT算法的调相信号Matlab实现

目录

一、什么是调相信号

二、调相信号的频谱是什么样的？

三、Matlab实现

---

上篇文章《全相位FFT算法Matlab实现》用的是一般的复数信号，但如果是调相信号，又该怎样用Matlab仿真实现呢？

---

一、什么是调相信号

看下面这张图。

---

二、调相信号的频谱是什么样的？

调相信号不同于一般的正弦信号，频谱分析需要用到贝塞尔函数进行展开，可参看如下论文，频谱见下图：《第一类贝塞尔函数在调频信号测量中的应用》

 

---

三、Matlab实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Zheng Wei, 2023/05/04
%%
%% 用途：如果信号频率f不等于fs/N的整数倍，FFT就会频谱泄露，计算的相位角就不对；
%%       全相位FFT法可以解决这个问题，但要注意，全相位FFT获取的是中心输入样点的相位值。比如
%%       当n=(-N+1:N-1)/fs时，apFFT所求的是n=0点的相位；
%%       当n=(1:2*N-1)/fs时， apFFT所求的是n=N点的相位。
%% 
%% 参考文献：《全相位FFT相位测量法》
%%
%% 全相位预处理过程：
%% 1、取一个N点窗函数w(n)；
%% 2、w(n) 自己对自己求卷积，得到2N-1的卷积窗；
%% 3、求2N-1的卷积窗的和；
%% 4、归一化卷积窗；
%% 5、将数据的1:2N-1项和归一化卷积窗相乘，得到加窗的2N-1项；
%% 6、将第1项和N+1项......第N-1项和第2N-1项相加，得到经过全相预处理N点序列；
%% 7、做N点FFT。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;        % 清理命令行
close all;  % 关闭所有图形窗口
clear all;  % 清理所有工作区变量
fclose all; % 关闭所有打开的文件

fc = 0;%0;% % 载波信号频率
fm = 55e6;%0;   % 调制信号频率
fd = 10e3;%0;   % 频率偏移
fs = 500e6;     % 采样率500Msps

N  = 1024;  % 如果f不等于fs/N的整数倍，就会造成频谱泄露，计算的相位角就不对

%% 当n =(-N+1:N-1)/fs时，apFFT所求中心点的相位：N_center = 0
%% 当n =(1:2*N-1)/fs时， apFFT所求中心点的相位：N_center = N
n =(-N+1:N-1)/fs;  % 文献上一般都是取这个范围
N_center = 0;
% n =(1:2*N-1)/fs;
% N_center = N;


%% 给个波形用于测试验证
r0 = pi/12;
k  = 2.0;%pi * 2.8 / 4; % 2.1991;
%x1_t = exp(j*(2*pi*fc*n+k*cos(2*pi*(fm-fd)*n+r0))+pi/20);
x1_t = exp(j*(2*pi*fc*n+k*cos(2*pi*(fm-fd)*n+r0)));
%x1_t   = exp(j*k*cos(2*pi*(fm-fd)*n+r0));


%% 手动计算出初始相位值，然后跟apFFT计算结果进行比对(可能会相差360°)
phase_x1 = r0 * 180/pi; % 角度
fprintf('手动计算的相位是：\t\t%d°!\n',phase_x1);


%% ==============================================================================
%% 以下是全相位FFT处理实现过程

%% 全相位预处理过程第1-4步：获取窗函数
win  = hanning(N)';
win1 = conv(win,win);
win2 = win1/sum(win1);


%% 全相位预处理过程第5步：将数据的1:2N-1项和归一化卷积窗相乘
Rx  = x1_t;
y1  = Rx .* win2;

%% 全相位预处理过程第6步：将第1项和N+1项......第N-1项和第2N-1项相加
y1a = y1(N:end) + [0 y1(1:N-1)];

%% 全相位预处理过程第7步：做N点FFT
Rx_fft = fft(y1a,N);

figure('Name','频谱');plot((0:N-1),abs(Rx_fft));title('全相位FFT之后的频谱');
xlim([-20,N+20]);

%% ==============================================================================
%% 全相位FFT处理完成之后，获取FFT最大峰值，其对应的相位值即是中心点相位值

Rx_fft = Rx_fft.* [ones(1,N/2) zeros(1,N/2)];

[fft_peak, fft_idx] = max(abs(Rx_fft));
max_Rx_fft     = Rx_fft(fft_idx);
apFFT_x1_rad   = angle(max_Rx_fft);    % 弧度
apFFT_x1_phase = apFFT_x1_rad * 180/pi; % 角度

fprintf('全相位FFT计算的相位是：  %d°!\n',apFFT_x1_phase);


%% ==============================================================================
%% 如下是直接使用FFT，可用于同apFFT处理效果进行对比

fft_x_n = fft(x1_t);
fft_x_n = fft_x_n.* [ones(1,N) zeros(1,N-1)];
[fft_x_peak, fft_x_idx] = max(abs(fft_x_n));
max_fft_x   = fft_x_n(fft_x_idx);
fft_x_rad   = angle(max_fft_x);    % 弧度
fft_x_phase = fft_x_rad * 180/pi   % 角度
figure('Name','频谱');plot((0:2*N-2),abs(fft_x_n));title('直接做FFT之后的频谱');
xlim([-20,2*N+20]);

两点特别说明：

1、特别需要注意调制度k值，要选取合适，我们选用了贝塞尔函数J1曲线；

2、FFT频谱只需要选择前面一半，否则结果可能相差180°；

3、计算结果相差90°，校正下即可。