如何用FFT提取非整数周期动态相位谱及200Hz相位参数

DevOps 小满

2026-7-19

问题描述

现有一频率为2kHz的光信号，经特定光路后生成包含两个不同拍频的光信号，该信号可简化为两个不同频率余弦信号的叠加——其中一个相位恒为0，另一个带有200Hz的正弦相位。我需要提取该200Hz正弦信号的相位值，即通过2kHz信号对200Hz相位进行采样，每周期提取一个相位值。但目前仅能使用每个周期的64%数据进行FFT分析，进而提取频谱相位，导致得到的相位结果存在误差。请问如何正确获取该200Hz相位的频率与幅值？

原测试代码

clc;clearvars;close all; 
v = 200;
Fs = 100e6; %采样率
fre = 2e3;%信号频率
count = 40;
n = 1.46; c = 3e8; w0 = 2*pi*3e8/1550e-9;
delta_v = 24e9; Tm =1/fre; points = Tm*Fs; t = 0:Tm/points:count*Tm-Tm/points; gamma =2*pi*delta_v/Tm; t=t';
duty = 0.8; ratio_start = 0.15; ratio_end = 0.95; %三角波占空比及截取范围 
len1 = 4; len2 = 7.8; len3 = 9; 
tau1 = 2*len1/c*n; tau2 = 2*len2/c*n; tau3 = 2*len3/c*n;

phase0 = 0; phase1 =0; phase2 =  3*sin(2*pi*v*t);  
E0 = (sawtooth(2*pi*fre*t,0.8)+1)/2;E1 = (sawtooth(2*pi*fre*t,0.8)+1)/2; E2 =(sawtooth(2*pi*fre*t,0.8)+1)/2; E3 = 0; 

E0t = E0.*exp(1i*(0.5*gamma*t.^2 + w0*t + phase0));
E1t = E1.*exp(1i*(0.5*gamma*(t-tau1).^2 + w0*(t-tau1) + phase1));
E2t = E2.*exp(1i*(0.5*gamma*(t-tau2).^2 + w0*(t-tau2) +phase1+ phase2));


Ut = (E0t + E1t).*conj(E0t + E1t) + (E0t + E2t).*conj(E0t + E2t) ;
Ut = Ut + wgn(length(E0t),1,-5) ;%添加噪声
Ut = Ut - mean(Ut);

eexx_valid_start = zeros(1,count-2);
eexx_valid_stop = zeros(1,count-2);
data = zeros(points*0.64,count-2);
x_t = zeros(points*0.64,count-2);
diff_phase = zeros(1,count-2);
locs = zeros(2,count-2);
start_point = 0;
N = 2^nextpow2(points*0.64);
f = Fs/N*(-N/2:1:N/2-1);f = f';
for i = 0:1:count-3
    eexx_valid_start(1,i+1) = floor(ratio_start*duty*points)+1+start_point+(i)*points;                        eexx_valid_stop(1,i+1) = floor(ratio_end*duty*points)+start_point+(i)*points;  
    data(:,i+1) = Ut(eexx_valid_start(1,i+1):eexx_valid_stop(1,i+1));
    x_t(:,i+1) = t(eexx_valid_start(1,i+1):eexx_valid_stop(1,i+1));
    F_data = fft(data(:,i+1),N)/(points*0.64)*2; F_data_shift = fftshift(F_data); 
    F_data_amp =  abs(F_data_shift);   
%     F_data_dB = 20*log10(F_data_amp);
    F_shift = F_data_shift;
    threshold = max(F_data_amp)/1000;
    F_shift(F_data_amp<threshold) = 0;
    phase = angle(F_shift);
    [piks,locs(:,i+1)] =       findpeaks(F_data_amp(N/2+100:N),f(N/2+100:N),'MinPeakHeight',0.6,'Annotate','extents');
    first_point = find(f == locs(1,i+1));   
    second_point = find(f == locs(2,i+1));
    diff_phase(1,i+1) = phase(second_point) - phase(first_point);
end
diff_phase_unwrap = unwrap(diff_phase(:)); 
%%
figure(1)
length2 = size(diff_phase_unwrap);
N2 = 16*2^nextpow2(length2(1));
F_y = fft(detrend(diff_phase_unwrap),N2)/length2(1)*2; F_y = fftshift(F_y); F_y = abs(F_y);
F_data_dB = 20*log10(F_y);
f2 = fre/N2*(-N2/2:1:N2/2-1);f2=f2'; 
plot(f2,F_data_dB);xlabel('频率/Hz'); ylabel('幅值/dB');xlim([0 2000]);

figure(2)
plot(detrend(diff_phase_unwrap))
xlabel('采样点'); ylabel('相位/rad');

解决思路与优化方案

1. 误差根源

当前用每个周期64%的截断数据做FFT，会引发频谱泄露，直接导致相位测量偏差；同时findpeaks定位频率点的方式易受噪声、旁瓣干扰，进一步放大相位误差。

2. 相位提取优化

(1) 加窗抑制频谱泄露

对截断数据先加汉宁窗（Hann）再做FFT，能有效降低旁瓣、减少泄露，同时修正幅值：

window = hann(length(data(:,i+1)));
data_windowed = data(:,i+1) .* window;
F_data = fft(data_windowed,N)/(sum(window)/2); % 加窗后幅值修正
F_data_shift = fftshift(F_data);

(2) 精确频率点插值

FFT的频率栅格是离散的，非整周期截断的真实频率不在栅格点上，用抛物线插值计算精确频率与对应相位：

% 假设找到的峰值索引为peak_idx
peak_idx = find(f == locs(1,i+1));
% 用相邻三点做抛物线插值求精确频率
amp_left = abs(F_data_shift(peak_idx-1));
amp_peak = abs(F_data_shift(peak_idx));
amp_right = abs(F_data_shift(peak_idx+1));
delta_f = (amp_left - amp_right)/(2*(2*amp_peak - amp_left - amp_right));
exact_f = f(peak_idx) + delta_f * Fs/N;
% 插值得到精确相位
exact_phase = angle(interp1(f, F_data_shift, exact_f));

(3) 固定相位参考

以恒相位的拍频为基准，每次计算时先将该拍频的相位对齐到0，再取另一个拍频的相位作为200Hz调制的采样值，确保相位差的参考一致性。

3. 200Hz信号参数提取

当得到准确的相位采样序列diff_phase_unwrap后：

频率：相位序列的采样率是2kHz（每2kHz周期采一个点），对去趋势后的序列做FFT，峰值对应的频率即为200Hz，需修正原代码中f2的计算逻辑：

Fs_phase = fre; % 相位序列采样率为2kHz
f2 = Fs_phase/N2*(-N2/2:1:N2/2-1);f2=f2';

幅值：相位序列为3*sin(2π*200*t)，FFT峰值的幅值需除以N/2（N为相位序列采样点数），即可得到实际幅值3。

内容的提问来源于stack exchange，提问作者future