Skip to content
Snippets Groups Projects
Commit 8b113a74 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Improved PT phase and number of bits allocation. 

Added phi_requant.
Added scale mode for satellite.
parent c0fefa56
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
applications/apertif/matlab/phase_tracking.m
+ 93
− 54
View file @ 8b113a74
...@@ -23,26 +23,32 @@ ...@@ -23,26 +23,32 @@
% Plot the phase tracking (PT) phase and phase erroras function of the % Plot the phase tracking (PT) phase and phase erroras function of the
% varying geometrical delay. % varying geometrical delay.
% Description: % Description:
% 1) Theory (see ASTRON-MEM-197)
% The geometrical delay tau_g varies due to the rotation of the earth. % The geometrical delay tau_g varies due to the rotation of the earth.
% The difference in geometrical delay between two antennas causes phase % The difference in geometrical delay between two antennas causes phase
% difference phi at the PT function: % difference phi at the PT function:
% %
% * phi_exact = w_RF * tau_g - w_n * tau_d % * 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 % The piecewise linear phase tracking tracks phi_exact using linear
% interpolation during each update interval T_update: % interpolation during each update interval T_update:
% %
% * phi_linear : piecewise linear intepolation double coefficients % * phi_linear : piecewise linear intepolation double coefficients
% * phi_quant : piecewise linear intepolation quantized 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 % 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 % scaling to determine the required number of bits for the PT
% coefficients. % coefficients.
% %
% By scaling the PT simulation can cover a half day, so theta 0-180 % 4) HDL reference data
% degrees. The scaling speeds up the earth rotation and therefore
% T_update needs to become smaller.
%
% Having a shorter T_update and larger T_sub is necessary to create % 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 % reference data for a HDL implementation of the PT. With these scaled
% parameters a HDL testbench that can still verify a few T_update % parameters a HDL testbench that can still verify a few T_update
...@@ -56,23 +62,28 @@ close all; ...@@ -56,23 +62,28 @@ close all;
% Control % Control
phi_error_max = 1; % choose maximum PT phase error in degrees 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 = 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 = -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 = 'no_scale '; % choose scaling of key geometrical delay parameters
%scale = false; scale = 'satellite';
scale = 'halfday ';
% Scaling % 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_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_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 scale_B = 10; % increase B to compensate for scale_f and preserve sufficient number of delay tracking DT steps Td
else elseif scale == 'satellite'
scale_w_e = 1; 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_f = 1; scale_T_update = 1; % do no compensate for increased scale_w_e, to check the impact on the PT accuracy
scale_B = 1;
end end
% Physical constants % Physical constants
...@@ -112,13 +123,13 @@ T_linear = sqrt(T_fine / w_e); ...@@ -112,13 +123,13 @@ T_linear = sqrt(T_fine / w_e);
% Time axis % Time axis
N_int = 800000; N_int = 800000;
T_int = N_int / f_sub; T_int = N_int / f_sub;
T_update = T_int / scale_w_e; % choose update interval equal to integration interval T_update = T_int / scale_T_update; % choose update interval equal to integration interval
nof_t_sub = round(N_int / scale_w_e); % number of subband time samples per update 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_i = (0:nof_t_sub-1);
t_u = t_i/f_sub; t_u = t_i/f_sub;
if T_linear < T_update 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 end
% Duration of the simulation % Duration of the simulation
...@@ -133,19 +144,22 @@ end ...@@ -133,19 +144,22 @@ end
T_sim_sub = nof_update * nof_t_sub; T_sim_sub = nof_update * nof_t_sub;
T_sim = nof_update * T_update; T_sim = nof_update * T_update;
% PT % PT phase and number of bits:
% * Number of bits: % <-- 360 degrees --->
% output phase phi : <-- phi_nof_bits --> % output phase phi : <-- phi_nof_bits -->
%
% input phase offset phi0 : <-- phi_nof_bits --><-- phi0_extra_bits --> % input phase offset phi0 : <-- phi_nof_bits --><-- phi0_extra_bits -->
% |----- N_phi ------||----- N_phi0_extra --|
% <-- phi0_nof_bits ------------------------> % <-- 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_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 % - 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 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; ...@@ -157,32 +171,44 @@ N_phi0_extra = 2^phi0_extra_bits;
N_phi0 = 2^phi0_nof_bits; N_phi0 = 2^phi0_nof_bits;
% - input phase step coefficient phi1 % - 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 d0_max = B / c;
N_phi1_fraction = 2^(phi0_extra_bits + phi1_fraction_nof_bits); d1_max = w_e * d0_max;
phi1_nof_bits = phi0_nof_bits + phi1_fraction_nof_bits; 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 = 2^phi1_nof_bits;
N_phi1_step = 2^phi1_step_nof_bits;
% Echo settings % Echo settings
sprintf (['Scale up factor for earth rotation rate = %d\n', ... sprintf (['Scale up factor for earth rotation rate = %d\n', ...
'Scale down factor for RF frequencies = %d\n', ... 'Scale down factor for update interval = %d\n', ...
'Scale up factor for baseline length = %d\n', ... 'Scale down factor for RF frequencies = %d\n', ...
'Scaled observation time for half day (12 hours) = %7.2f [s]\n', ... 'Scale up factor for baseline length = %d\n', ...
'Simulated observation time = %7.2f [s]\n', ... 'Scaled observation time for half day (12 hours) = %7.2f [s]\n', ...
'Subband RF frequency = %7.2f [MHz]\n', ... 'Simulated observation time = %7.2f [s]\n', ...
'DT delay step Td = %7.3f [ns]\n', ... 'Subband RF frequency = %7.2f [MHz]\n', ...
'Maximum nof delay steps Td = %d\n', ... 'DT delay step Td = %7.3f [ns]\n', ...
'Maximum fringe frequency = %7.2f [Hz]\n', ... 'Maximum nof delay steps Td = %d\n', ...
'Maximum linear update interval = %7.5f [s]\n', ... 'Maximum fringe frequency = %7.2f [Hz]\n', ...
'Linear update interval = %7.5f [s]\n', ... 'Maximum linear update interval = %7.5f [s]\n', ...
'Nof update intervals = %d\n', ... 'Linear update interval = %7.5f [s]\n', ...
'Nof subband time samples per update interval = %d\n', ... 'Nof update intervals = %d\n', ...
'Nof subband time samples in total simulation = %d\n', ... 'Nof subband time samples per update interval = %d\n', ...
'Maximum PT phase error = %7.3f [degrees]\n', ... 'Nof subband time samples in total simulation = %d\n', ...
'Nof bits for 360 degrees output phase = %d\n', ... 'Maximum phase step phi1 = %f [degrees]\n', ...
'Nof bits for piecewise linear input phase offset = %d\n', ... 'Maximum PT phase error = %7.3f [degrees]\n', ...
'Nof bits for piecewise linear input phase step = %d\n'], ... 'Nof bits for 360 degrees input phase phi0 = %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,... 'Nof bits for input interpolation interval = %d\n', ...
phi_error_max, phi_nof_bits, phi0_nof_bits, phi1_nof_bits) '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) % Simulate the phase tracking (PT)
r_theta0 = g_theta0 * pi/180; r_theta0 = g_theta0 * pi/180;
...@@ -190,6 +216,7 @@ DT = []; ...@@ -190,6 +216,7 @@ DT = [];
phi_exact = []; phi_exact = [];
phi_linear = []; phi_linear = [];
phi_quant = []; phi_quant = [];
phi_requant = [];
hdl_phi_coefficients = []; hdl_phi_coefficients = [];
hdl_phi_phases = []; hdl_phi_phases = [];
for up = 0:nof_update-1 for up = 0:nof_update-1
...@@ -220,16 +247,18 @@ for up = 0:nof_update-1 ...@@ -220,16 +247,18 @@ for up = 0:nof_update-1
phi1 = g_RF * d1 * T_sub; phi1 = g_RF * d1 * T_sub;
phi1_q = round(N_phi1 * phi1/360); phi1_q = round(N_phi1 * phi1/360);
% - requantize the PT output % - requantize the PT output
phi_q = round(phi0_q / N_phi0_extra + phi1_q * t_i / N_phi1_fraction); 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]; 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 % Keep reference data for HDL simulation
% - input coefficients per update interval % - input coefficients per update interval
hdl_phi0_q = rem(phi0_q, N_phi0); % effectively store the phase offset coefficient modulo 2pi 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]]; hdl_phi_coefficients = [hdl_phi_coefficients; [hdl_phi0_q hdl_phi1_q]];
% - expected output phases per update interval % - expected output phases per update interval
hdl_phi_phases = [hdl_phi_phases; phi_q]; hdl_phi_phases = [hdl_phi_phases; phi_rq];
end end
g_theta_end = (r_theta0 + w_e * T_update * nof_update) * 180/pi; 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_ ...@@ -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]']); xlabel(['Update intervals, T update = ' num2str(T_update) ' [s]']);
ylabel(['Phase error [degrees]']); ylabel(['Phase error [degrees]']);
grid on; 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;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment