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

Added quantization, proper scaling and many small improvements.

parent 84692a42
Branches
No related tags found
No related merge requests found
...@@ -21,8 +21,27 @@ ...@@ -21,8 +21,27 @@
% Author: W. Lubberhuizen, 2015 (top.m for LOFAR Station, uses classes) % Author: W. Lubberhuizen, 2015 (top.m for LOFAR Station, uses classes)
% Eric Kooistra, 2016 % Eric Kooistra, 2016
% Purpose : Model data path for delay tracking wiht polyphase filterbank % Purpose : Model data path for delay tracking with polyphase filterbank
% Description : % Description :
% * General
% The data path (DP) is modelled per block of data. The block of data
% are counted by the block sequence number that thus acts as a
% timestamp. At each BSN there is the possibility to apply monitoring
% and control (MAC). The data path consist of DSP functions. Each DSP
% function expects normalized full scale input (range -1:1) and it also
% normalized full scale output. Dependend on the tb.model the DSP
% functions outputs floating point or fixed point data that is quantized
% to the specified number of bits. The output full scale range remains
% -1:1.
% * Fringe stopping
% The delay step of the fringe stopping causes a phase step in the
% subband phase. This phase step occurs after nof_taps/2. If nof_taps is
% odd then the phase step occurs within one subband periods, else when
% even it takes two subband periods. Hence it may be appropriate to
% delay the application of the new PT coefficients by floor(nof_taps/2)
% subband periods after the delay step of the delay tracking. Still the
% delay step also causes amplitude distortion and aliasing into all
% other subbands, which cannot be mitigated by delaying the phase step.
clear all; clear all;
close all; close all;
...@@ -31,19 +50,42 @@ fig=0; ...@@ -31,19 +50,42 @@ fig=0;
tb.nof_complex = 2; tb.nof_complex = 2;
tb.nof_subbands = 512; tb.nof_subbands = 512;
tb.bsn_init = 0; tb.bsn_init = 0;
tb.subband_nr = 101; % subband range 0:511, can be fraction tb.subband_nr = 113; % subband range 0:511, can be fraction
%tb.subband_nr = 113.3;
tb.subband_fft_size = tb.nof_complex*tb.nof_subbands; % subband filterbank real FFT tb.subband_fft_size = tb.nof_complex*tb.nof_subbands; % subband filterbank real FFT
tb.block_size = tb.subband_fft_size; % number of samples per block tb.nof_blocks = 61;
tb.nof_blocks = 53; tb.model = 'floating point';
tb.model = 'fixed point';
% DP quantization, width 0 is use double, width > 0 is use nof bits for data
if strcmp(tb.model, 'floating point')
ctrl_wg.data_w = 0;
ctrl_pfir_subband.coeff_w = 0;
ctrl_pfir_subband.scale_w = 0;
ctrl_pfir_subband.data_w = 0;
ctrl_pfft_subband.data_w = 0;
ctrl_spec_subband.db_low = -150;
else
ctrl_wg.data_w = 8;
ctrl_pfir_subband.coeff_w = 9; % bits/coefficient = 1 sign + (w-1) mantissa
ctrl_pfir_subband.scale_w = 1; % scale data by 2^scale before quantization to account for PFIR overshoot
ctrl_pfir_subband.data_w = 16;
ctrl_pfft_subband.data_w = 16;
ctrl_spec_subband.data_w = 10;
ctrl_spec_subband.db_low = floor(-20 - 6.02 * ctrl_spec_subband.data_w);
end
% Waveform generator % Waveform generator
ctrl_wg.block_size = tb.block_size; ctrl_wg.block_size = tb.subband_fft_size;
ctrl_wg.ampl = 1; ctrl_wg.ampl = 1; % full scale is 1, ampl > 1 causes analogue clipping
ctrl_wg.offset = 0; % DC offset
ctrl_wg.freq = tb.subband_nr/tb.subband_fft_size; % normalized fs ctrl_wg.freq = tb.subband_nr/tb.subband_fft_size; % normalized fs
ctrl_wg.phase = 0; % normalized 2pi ctrl_wg.phase = 0; % normalized 2pi
% Delay tracking % Delay tracking
ctrl_dt.block_size = tb.block_size; ctrl_dt.block_size = tb.subband_fft_size;
ctrl_dt.buffer = zeros(1, 2*ctrl_dt.block_size); ctrl_dt.buffer = zeros(1, 2*ctrl_dt.block_size);
ctrl_dt.step = 4; % Delay step is Psub = 4 factor, for 4 samples per clock ctrl_dt.step = 4; % Delay step is Psub = 4 factor, for 4 samples per clock
ctrl_dt.dt = 0; % Delay setting in +- number of ctrl_dt.step time samples of ADC or WG ctrl_dt.dt = 0; % Delay setting in +- number of ctrl_dt.step time samples of ADC or WG
...@@ -59,15 +101,15 @@ ctrl_pfir_subband.r_pass = 0.001; ...@@ -59,15 +101,15 @@ ctrl_pfir_subband.r_pass = 0.001;
ctrl_pfir_subband.r_stop = 0.0001; ctrl_pfir_subband.r_stop = 0.0001;
ctrl_pfir_subband.hp_factor = 1.050; % Adjust channel half power bandwidth ctrl_pfir_subband.hp_factor = 1.050; % Adjust channel half power bandwidth
ctrl_pfir_subband.hp_factor = 1; ctrl_pfir_subband.hp_factor = 1;
ctrl_pfir_subband.Wchan = ctrl_pfir_subband.hp_factor / tb.subband_fft_size; % Channel bandwidth ctrl_pfir_subband.BWchan = ctrl_pfir_subband.hp_factor / tb.subband_fft_size; % Channel bandwidth
ctrl_pfir_subband.coeff_w = 0; % bits/coefficient = 1 sign + (W-1) mantissa, 0 for floating point
ctrl_pfir_subband.config.design = 'fir1'; % 'fir1', 'fircls1', 'lofar_file' ctrl_pfir_subband.config.design = 'fir1'; % 'fir1', 'fircls1', 'lofar_file'
ctrl_pfir_subband.config.design = 'fircls1';
ctrl_pfir_subband.config.design_flag = 'trace'; % 'trace' ctrl_pfir_subband.config.design_flag = 'trace'; % 'trace'
ctrl_pfir_subband.config.interpolate = 'interpft'; % 'resample', 'fourier', 'interpft' ctrl_pfir_subband.config.interpolate = 'interpft'; % 'resample', 'fourier', 'interpft'
coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,.... coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,....
ctrl_pfir_subband.nof_taps, ... ctrl_pfir_subband.nof_taps, ...
ctrl_pfir_subband.Wchan, ... ctrl_pfir_subband.BWchan, ...
ctrl_pfir_subband.r_pass, ... ctrl_pfir_subband.r_pass, ...
ctrl_pfir_subband.r_stop, ... ctrl_pfir_subband.r_stop, ...
ctrl_pfir_subband.coeff_w, ... ctrl_pfir_subband.coeff_w, ...
...@@ -77,19 +119,21 @@ coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,.... ...@@ -77,19 +119,21 @@ coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,....
%ctrl_pfir_subband.coeff = flipud(coeff); %ctrl_pfir_subband.coeff = flipud(coeff);
ctrl_pfir_subband.coeff = coeff; ctrl_pfir_subband.coeff = coeff;
ctrl_pfir_subband.Zdelays = zeros(ctrl_pfir_subband.downsample_factor, ctrl_pfir_subband.nof_taps-1); ctrl_pfir_subband.Zdelays = zeros(ctrl_pfir_subband.downsample_factor, ctrl_pfir_subband.nof_taps-1);
ctrl_pfir_subband.gain = sum(ctrl_pfir_subband.coeff(:)) / ctrl_pfir_subband.downsample_factor;
% Subband FFT parameters % Subband FFT parameters
ctrl_pfft_subband.fft_size = tb.subband_fft_size; ctrl_pfft_subband.fft_size = tb.subband_fft_size;
ctrl_pfft_subband.complex = false; ctrl_pfft_subband.complex = false;
ctrl_pfft_subband.gain = ctrl_pfft_subband.fft_size;
% Run the data path processing % Run the data path processing (one block per row)
t_start = cputime; t_start = cputime;
data_dt = []; data_dt = zeros(tb.nof_blocks, tb.subband_fft_size);
data_pfir_subband = []; data_pfir_subband = zeros(tb.nof_blocks, tb.subband_fft_size);
data_pfft_subband = []; data_pfft_subband = zeros(tb.nof_blocks, tb.nof_subbands);
for block_nr = 1:tb.nof_blocks for bi = 1:tb.nof_blocks
% Control % Control
if ctrl_bsn_source.bsn == 32 if ctrl_bsn_source.bsn == 30
ctrl_dt.dt = ctrl_dt.dt - ctrl_dt.step; % Apply delay step ctrl_dt.dt = ctrl_dt.dt - ctrl_dt.step; % Apply delay step
end end
...@@ -101,15 +145,15 @@ for block_nr = 1:tb.nof_blocks ...@@ -101,15 +145,15 @@ for block_nr = 1:tb.nof_blocks
[ block_pfft_subband] = pfft( ctrl_pfft_subband, block_pfir_subband); [ block_pfft_subband] = pfft( ctrl_pfft_subband, block_pfir_subband);
% Capture data at each DP interface % Capture data at each DP interface
data_dt = [data_dt block_dt]; data_dt(bi, :) = block_dt;
data_pfir_subband = [data_pfir_subband block_pfir_subband]; data_pfir_subband(bi, :) = block_pfir_subband;
data_pfft_subband = [data_pfft_subband block_pfft_subband]; data_pfft_subband(bi, :) = block_pfft_subband;
end end
t_stop = cputime; t_stop = cputime;
disp(sprintf('Total processing time: %f seconds', t_stop-t_start)); disp(sprintf('Total processing time: %f seconds', t_stop-t_start));
% Plot data path results % Plot data path results
ts = (0:tb.block_size*tb.nof_blocks-1)/tb.block_size; % time of ADC / WG samples in block periods ts = (0:tb.subband_fft_size*tb.nof_blocks-1)/tb.subband_fft_size; % time of ADC / WG samples in block periods
tsub = (0:tb.nof_subbands*tb.nof_blocks-1)/tb.nof_subbands; % time in subband periods tsub = (0:tb.nof_subbands*tb.nof_blocks-1)/tb.nof_subbands; % time in subband periods
tblock = (0:tb.nof_blocks-1); % time in blocks tblock = (0:tb.nof_blocks-1); % time in blocks
...@@ -117,53 +161,79 @@ tblock = (0:tb.nof_blocks-1); % time in blocks ...@@ -117,53 +161,79 @@ tblock = (0:tb.nof_blocks-1); % time in blocks
fig=fig+1; fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]); figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig); figure(fig);
plot(ts, data_dt) data = data_dt';
ylim(1.2*[-ctrl_wg.ampl ctrl_wg.ampl]); plot(ts, data(:))
title('Delay tracking output data'); ylim([-1.3 1.3]);
title(sprintf('Delay tracking output data (WG subband %6.3f)', tb.subband_nr));
xlabel('Time [T bsn]'); xlabel('Time [T bsn]');
ylabel('Voltage'); ylabel('Voltage');
grid on; grid on;
% Plot FIR output % Plot PFIR transfer function
h = ctrl_pfir_subband.coeff;
h = h(:);
NL = ctrl_pfir_subband.nof_coefficients;
L = ctrl_pfir_subband.nof_taps;
hf_abs = abs(fftshift(fft(h / sum(h), NL)));
hi = 1:NL; % coefficients index
fx = (hi - NL/2-1) / L; % frequency axis in subband units
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(fx, db(hf_abs)); % db() = 20*log10 for voltage
%xlim([-3 3]);
grid on;
title(['FIR filter transfer function for ', ctrl_pfir_subband.config.design]);
xlabel('Frequency [channels]');
ylabel('Magnitude [dB]');
% Plot PFIR output
fig=fig+1; fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]); figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig); figure(fig);
plot(ts, data_pfir_subband) data = data_pfir_subband';
title('Subband polyphase FIR filter output - FFT input data'); plot(ts, data(:))
title(sprintf('Subband polyphase FIR filter output - FFT input data (WG subband %6.3f)', tb.subband_nr));
ylim([-1.3 1.3]);
xlabel('Time [T bsn]'); xlabel('Time [T bsn]');
ylabel('Voltage'); ylabel('Voltage');
grid on; grid on;
% Plot FFT subbands spectrum and phase % Plot PFFT subbands spectrum and phase
ampl = abs(data_pfft_subband);
phase = angle(data_pfft_subband)*180/pi;
noise = ampl < 0.1*max(ampl(:));
phase(noise) = 0; % force phase of too small signals to 0
fig=fig+1; fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]); figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig); figure(fig);
subplot(2,1,1); subplot(2,1,1);
plot(tsub, abs(data_pfft_subband)) data = ampl';
plot(tsub, data(:))
title('Subband data - amplitude'); title('Subband data - amplitude');
xlabel('Time [T bsn]'); xlabel('Time [T bsn]');
ylabel('Voltage'); ylabel('Voltage');
grid on; grid on;
subplot(2,1,2); subplot(2,1,2);
% force phase of too small signals to 0 data = phase';
ampl = abs(data_pfft_subband); plot(tsub, data(:))
noise = ampl<0.1*max(ampl); ylim([-180 180])
phase = angle(data_pfft_subband)*180/pi;
phase(noise) = 0;
plot(tsub, phase)
title('Subband data - phase'); title('Subband data - phase');
xlabel('Time [T bsn]'); xlabel('Time [T bsn]');
ylabel('Phase [degrees]'); ylabel('Phase [degrees]');
grid on; grid on;
% Plot phase step due to DT step for the subband that is set in the WG % Plot phase step due to DT step for the subband that is set in the WG
nr = round(tb.subband_nr);
phase = angle(data_pfft_subband(:, nr+1))*180/pi;
fig=fig+1; fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]); figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig); figure(fig);
nr = round(tb.subband_nr);
phase = angle(downsample(data_pfft_subband, tb.nof_subbands, nr))*180/pi;
plot(tblock, phase) plot(tblock, phase)
title(sprintf('Subband phase for subband %f in range 0:%d', nr, tb.nof_subbands-1)); ylim([-180 180])
title(sprintf('Subband phase for subband %d in range 0:%d', nr, tb.nof_subbands-1));
xlabel('Time [T bsn]'); xlabel('Time [T bsn]');
ylabel('Phase [degrees]'); ylabel('Phase [degrees]');
grid on; grid on;
...@@ -173,18 +243,17 @@ fig=fig+1; ...@@ -173,18 +243,17 @@ fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]); figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig); figure(fig);
data = abs(data_pfft_subband); data = abs(data_pfft_subband);
data = data / max(data); data_abs_max = max(data(:));
data = db(data); data = db(data); % no need to scale data, range is already normalized
x = floor(min(data(data~=-Inf))); x = ctrl_spec_subband.db_low;
%x = -120; %x = floor(min(data(data~=-Inf)))
%data(data<x) = x; data(data<x) = x;
max(data)
min(data)
data = reshape(data, tb.nof_subbands, tb.nof_blocks);
color_range = size(colormap, 1);
mymap = jet(-x); mymap = jet(-x);
%mymap(1,:) = [0 0 0]; %mymap(1,:) = [0 0 0]; % force black for dB(0)
colormap(mymap); colormap(mymap);
imagesc(data); imagesc(data',[x 0]);
colorbar; colorbar;
title(sprintf('Subband spectogram (max value = %f = %f dB)', data_abs_max, db(data_abs_max)));
xlabel('Time [T bsn]');
ylabel('Subband');
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment