目录
上篇文章《全相位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°,校正下即可。