diff --git a/applications/apertif/matlab/phase_tracking.m b/applications/apertif/matlab/phase_tracking.m
index 09ba585fbd7b61f2edddce27216b4ccd0cccf957..c889ef3a33dbbdaa3ffbb93360c1ad23cd14d881 100644
--- a/applications/apertif/matlab/phase_tracking.m
+++ b/applications/apertif/matlab/phase_tracking.m
@@ -23,26 +23,32 @@
% Plot the phase tracking (PT) phase and phase erroras function of the
% varying geometrical delay.
% Description:
+% 1) Theory (see ASTRON-MEM-197)
% The geometrical delay tau_g varies due to the rotation of the earth.
% The difference in geometrical delay between two antennas causes phase
% difference phi at the PT function:
%
% * phi_exact = w_RF * tau_g - w_n * tau_d
%
+% By scaling the PT simulation can cover a half day, so theta 0-180
+% degrees. The scaling speeds up the earth rotation and therefore
+% T_update needs to become smaller.
+%
+% 2) Piecewise linear approximation
% The piecewise linear phase tracking tracks phi_exact using linear
% interpolation during each update interval T_update:
%
-% * phi_linear : piecewise linear intepolation double coefficients
-% * phi_quant : piecewise linear intepolation quantized coefficients
+% * phi_linear : piecewise linear intepolation double coefficients
+%
+% 3) Quantization
+% * phi_quant : piecewise linear intepolation quantized coefficients
+% * phi_requant : piecewise linear intepolation requantized coefficients
%
% The required accuracy for the PT is set via phi_error_max. Use no
% scaling to determine the required number of bits for the PT
% coefficients.
%
-% By scaling the PT simulation can cover a half day, so theta 0-180
-% degrees. The scaling speeds up the earth rotation and therefore
-% T_update needs to become smaller.
-%
+% 4) HDL reference data
% Having a shorter T_update and larger T_sub is necessary to create
% reference data for a HDL implementation of the PT. With these scaled
% parameters a HDL testbench that can still verify a few T_update
@@ -56,23 +62,28 @@ close all;
% Control
phi_error_max = 1; % choose maximum PT phase error in degrees
phi0_extra_bits = 0; % no extra bits are needed to achieve phi_error_max, use e.g. 8 extra bits to have phi_quant approach phi_linear
-%phi0_extra_bits = 8;
+phi0_extra_bits = 7;
g_theta0 = -90; % choose start observation angle in degrees (-90 is east, 0 is up, +90 is west)
-g_theta0 = 0;
+%g_theta0 = 0;
-scale = true; % choose scaling of key geometrical delay parameters or not
-%scale = false;
+scale = 'no_scale '; % choose scaling of key geometrical delay parameters
+scale = 'satellite';
+scale = 'halfday ';
% Scaling
-if scale
+scale_w_e = 1; % default no scaling
+scale_T_update = 1;
+scale_f = 1;
+scale_B = 1;
+if scale == 'halfday '
scale_w_e = 10000; % speed up earth rotation to effectively make shorter day
+ scale_T_update = scale_w_e;
scale_f = 10; % decrease f_RF to have larger T_sub unit time steps
scale_B = 10; % increase B to compensate for scale_f and preserve sufficient number of delay tracking DT steps Td
-else
- scale_w_e = 1;
- scale_f = 1;
- scale_B = 1;
+elseif scale == 'satellite'
+ scale_w_e = 15; % speed up earth rotation to effectively model the angular speed of a satellite (e.g. 24 h / 15 = 1.6 hour orbit)
+ scale_T_update = 1; % do no compensate for increased scale_w_e, to check the impact on the PT accuracy
end
% Physical constants
@@ -112,13 +123,13 @@ T_linear = sqrt(T_fine / w_e);
% Time axis
N_int = 800000;
T_int = N_int / f_sub;
-T_update = T_int / scale_w_e; % choose update interval equal to integration interval
-nof_t_sub = round(N_int / scale_w_e); % number of subband time samples per update interval
+T_update = T_int / scale_T_update; % choose update interval equal to integration interval
+nof_t_sub = round(N_int / scale_T_update); % number of subband time samples per update interval
t_i = (0:nof_t_sub-1);
t_u = t_i/f_sub;
if T_linear < T_update
- sprintf('Error: the update interval is too long for piecewise linear interpolation !!!\n')
+ sprintf('Warning: the update interval is too long for piecewise linear interpolation !!!\n')
end
% Duration of the simulation
@@ -133,19 +144,22 @@ end
T_sim_sub = nof_update * nof_t_sub;
T_sim = nof_update * T_update;
-% PT
-% * Number of bits:
+% PT phase and number of bits:
+% <-- 360 degrees --->
% output phase phi : <-- phi_nof_bits -->
+%
% input phase offset phi0 : <-- phi_nof_bits --><-- phi0_extra_bits -->
+% |----- N_phi ------||----- N_phi0_extra --|
% <-- phi0_nof_bits ------------------------>
-% input phase step phi1 : <-- phi0_nof_bits ------------------------><-- phi1_fraction_nof_bits -->
-% <-- phi1_nof_bits ------------------------------------------------------>
-% * Corresponding ranges:
-% |----- N_phi ------|
-% |----- N_phi0_extra ---|
% |----- N_phi0 ----------------------------|
-% |----- N_phi1_fraction ------------------------------|
-% |----- N_phi1 ----------------------------------------------------------|
+%
+% input phase step phi1 : <-- phi_nof_bits --><-- phi0_extra_bits --><-- phi1_interpolate_nof_bits -->
+% <-- phi1_fraction_nof_bits ---------------------------->
+% |----- N_phi1_fraction --------------------------------|
+% <-- phi1_nof_bits --------------------------------------------------------->
+% |----- N_phi1 -------------------------------------------------------------|
+% <-- phi1_backoff_nof_bits --><-- phi1_step_nof_bits ----------------------->
+% |----- N_phi1_step ---------------------------|
% - output phase phi
phi_nof_bits = ceil(log2(360/phi_error_max)); % nof bits to represent 0:360, used for phi -> (Re,Im) look up table address range in HDL
@@ -157,32 +171,44 @@ N_phi0_extra = 2^phi0_extra_bits;
N_phi0 = 2^phi0_nof_bits;
% - input phase step coefficient phi1
-phi1_fraction_nof_bits = ceil(log2(nof_t_sub/phi_error_max)); % nof bits to preserve phi_error_max after adding nof_t_sub
-N_phi1_fraction = 2^(phi0_extra_bits + phi1_fraction_nof_bits);
-phi1_nof_bits = phi0_nof_bits + phi1_fraction_nof_bits;
+d0_max = B / c;
+d1_max = w_e * d0_max;
+phi1_max = g_RF * d1_max * T_sub;
+phi1_interpolate_nof_bits = ceil(log2(nof_t_sub/phi_error_max)); % nof bits to preserve phi_error_max after adding nof_t_sub
+phi1_fraction_nof_bits = phi0_extra_bits + phi1_interpolate_nof_bits; % nof bits to preserve phi_error_max after adding nof_t_sub and phi0_extra_bits
+phi1_nof_bits = phi_nof_bits + phi1_fraction_nof_bits;
+phi1_backoff_nof_bits = floor(log2(360 / phi1_max)) - 1; % -1 for sign bit of phase step phi1
+phi1_step_nof_bits = phi1_nof_bits - phi1_backoff_nof_bits;
+N_phi1_fraction = 2^phi1_fraction_nof_bits;
N_phi1 = 2^phi1_nof_bits;
+N_phi1_step = 2^phi1_step_nof_bits;
% Echo settings
-sprintf (['Scale up factor for earth rotation rate = %d\n', ...
- 'Scale down factor for RF frequencies = %d\n', ...
- 'Scale up factor for baseline length = %d\n', ...
- 'Scaled observation time for half day (12 hours) = %7.2f [s]\n', ...
- 'Simulated observation time = %7.2f [s]\n', ...
- 'Subband RF frequency = %7.2f [MHz]\n', ...
- 'DT delay step Td = %7.3f [ns]\n', ...
- 'Maximum nof delay steps Td = %d\n', ...
- 'Maximum fringe frequency = %7.2f [Hz]\n', ...
- 'Maximum linear update interval = %7.5f [s]\n', ...
- 'Linear update interval = %7.5f [s]\n', ...
- 'Nof update intervals = %d\n', ...
- 'Nof subband time samples per update interval = %d\n', ...
- 'Nof subband time samples in total simulation = %d\n', ...
- 'Maximum PT phase error = %7.3f [degrees]\n', ...
- 'Nof bits for 360 degrees output phase = %d\n', ...
- 'Nof bits for piecewise linear input phase offset = %d\n', ...
- 'Nof bits for piecewise linear input phase step = %d\n'], ...
- scale_w_e, scale_f, scale_B, T_halfday, T_sim, f_RF/1e6 , Td*1e9, k_max, f_fringe_max , T_linear, T_update, nof_update, nof_t_sub, T_sim_sub,...
- phi_error_max, phi_nof_bits, phi0_nof_bits, phi1_nof_bits)
+sprintf (['Scale up factor for earth rotation rate = %d\n', ...
+ 'Scale down factor for update interval = %d\n', ...
+ 'Scale down factor for RF frequencies = %d\n', ...
+ 'Scale up factor for baseline length = %d\n', ...
+ 'Scaled observation time for half day (12 hours) = %7.2f [s]\n', ...
+ 'Simulated observation time = %7.2f [s]\n', ...
+ 'Subband RF frequency = %7.2f [MHz]\n', ...
+ 'DT delay step Td = %7.3f [ns]\n', ...
+ 'Maximum nof delay steps Td = %d\n', ...
+ 'Maximum fringe frequency = %7.2f [Hz]\n', ...
+ 'Maximum linear update interval = %7.5f [s]\n', ...
+ 'Linear update interval = %7.5f [s]\n', ...
+ 'Nof update intervals = %d\n', ...
+ 'Nof subband time samples per update interval = %d\n', ...
+ 'Nof subband time samples in total simulation = %d\n', ...
+ 'Maximum phase step phi1 = %f [degrees]\n', ...
+ 'Maximum PT phase error = %7.3f [degrees]\n', ...
+ 'Nof bits for 360 degrees input phase phi0 = %d\n', ...
+ 'Nof bits for input interpolation interval = %d\n', ...
+ 'Nof bits for input phase step phi1 = %d\n', ...
+ 'Nof bits for 360 degrees internal phase = %d\n', ...
+ 'Nof bits for 360 degrees internal phase fraction = %d\n', ...
+ 'Nof bits for 360 degrees output phase phi = %d\n'], ...
+ scale_w_e, scale_T_update, scale_f, scale_B, T_halfday, T_sim, f_RF/1e6 , Td*1e9, k_max, f_fringe_max , T_linear, T_update, nof_update, nof_t_sub, T_sim_sub,...
+ phi1_max, phi_error_max, phi0_nof_bits, phi1_interpolate_nof_bits, phi1_step_nof_bits, phi1_nof_bits, phi1_fraction_nof_bits, phi_nof_bits)
% Simulate the phase tracking (PT)
r_theta0 = g_theta0 * pi/180;
@@ -190,6 +216,7 @@ DT = [];
phi_exact = [];
phi_linear = [];
phi_quant = [];
+phi_requant = [];
hdl_phi_coefficients = [];
hdl_phi_phases = [];
for up = 0:nof_update-1
@@ -220,16 +247,18 @@ for up = 0:nof_update-1
phi1 = g_RF * d1 * T_sub;
phi1_q = round(N_phi1 * phi1/360);
% - requantize the PT output
- phi_q = round(phi0_q / N_phi0_extra + phi1_q * t_i / N_phi1_fraction);
- phi_quant = [phi_quant phi_q * 360/N_phi];
+ phi_q = phi0_q / N_phi0_extra + phi1_q * t_i / N_phi1_fraction; % quantized real output phi_q
+ phi_quant = [phi_quant phi_q * 360/N_phi]; % use real phi_q output to see effect of increasing phi0_extra_bits
+ phi_rq = round(phi_q); % requantized integer output phi_rq
+ phi_requant = [phi_requant phi_rq * 360/N_phi]; % use integer phi_rq output to see actual requantized output
% Keep reference data for HDL simulation
% - input coefficients per update interval
hdl_phi0_q = rem(phi0_q, N_phi0); % effectively store the phase offset coefficient modulo 2pi
- hdl_phi1_q = rem(phi1_q, N_phi1); % effectively store the phase step coefficient modulo 2pi, no need to do modulo N_phi1,use rem to preserve sign
+ hdl_phi1_q = rem(phi1_q, N_phi1_step); % effectively store the phase step coefficient modulo maximum phase step, use rem to preserve sign
hdl_phi_coefficients = [hdl_phi_coefficients; [hdl_phi0_q hdl_phi1_q]];
% - expected output phases per update interval
- hdl_phi_phases = [hdl_phi_phases; phi_q];
+ hdl_phi_phases = [hdl_phi_phases; phi_rq];
end
g_theta_end = (r_theta0 + w_e * T_update * nof_update) * 180/pi;
@@ -297,3 +326,13 @@ title(['Phi exact - phi quant, for theta = ' num2str(g_theta0) ' to ' num2str(g_
xlabel(['Update intervals, T update = ' num2str(T_update) ' [s]']);
ylabel(['Phase error [degrees]']);
grid on;
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(tPhi, phi_exact-phi_requant)
+%axis([0 nof_update -phi_error_max phi_error_max]);
+title(['Phi exact - phi requant, for theta = ' num2str(g_theta0) ' to ' num2str(g_theta_end) ' [degrees]']);
+xlabel(['Update intervals, T update = ' num2str(T_update) ' [s]']);
+ylabel(['Phase error [degrees]']);
+grid on;