How to extract the dynamic phase spectrum of non-integer cycles with FFT

Question

I have an optical signal with a frequency of 2K, and then after a certain optical path, an optical signal consisting of two different beat frequencies is generated. This optical signal can be simply seen as a superposition of two cosine signals of different frequencies, one of which has a constant phase of 0 and the other a sine signal with a phase of 200Hz.

I need to extract the phase value of this 200Hz sine, i.e., by sampling the 200Hz phase with a 2K signal, and extracting one phase value per cycle. But unfortunately, I can only use 64% of each cycle at a time, perform an FFT on it, and then extract the phase of the spectrum. This makes the phase results I get incorrect. How can I get the frequency and amplitude of that 200Hz phase correctly.

clc;clearvars;close all; 
v = 200;
Fs = 100e6; %Sampling rate
fre = 2e3;%singal frequency
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; %Duty cycle of triangular waveform and interception range 
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) ;%noise
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('f/Hz'); ylabel('Amp/dB');xlim([0 2000]);

figure(2)
plot(detrend(diff_phase_unwrap))
xlabel('f/Hz'); ylabel('phase/rad');

enter image description hereenter image description here