diff --git a/applications/apertif/matlab/fringe_stopping.m b/applications/apertif/matlab/fringe_stopping.m index 920e1315c8a74f2cab02d5db790e7b9def98d366..25694a68c444362865ee3c4ef05cd70c73904dc3 100644 --- a/applications/apertif/matlab/fringe_stopping.m +++ b/applications/apertif/matlab/fringe_stopping.m @@ -21,8 +21,27 @@ % Author: W. Lubberhuizen, 2015 (top.m for LOFAR Station, uses classes) % Eric Kooistra, 2016 -% Purpose : Model data path for delay tracking wiht polyphase filterbank +% Purpose : Model data path for delay tracking with polyphase filterbank % 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; close all; @@ -31,19 +50,42 @@ fig=0; tb.nof_complex = 2; tb.nof_subbands = 512; 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.block_size = tb.subband_fft_size; % number of samples per block -tb.nof_blocks = 53; +tb.nof_blocks = 61; +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 -ctrl_wg.block_size = tb.block_size; -ctrl_wg.ampl = 1; +ctrl_wg.block_size = tb.subband_fft_size; +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.phase = 0; % normalized 2pi % 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.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 @@ -59,15 +101,15 @@ ctrl_pfir_subband.r_pass = 0.001; 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; -ctrl_pfir_subband.Wchan = 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.BWchan = ctrl_pfir_subband.hp_factor / tb.subband_fft_size; % Channel bandwidth 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.interpolate = 'interpft'; % 'resample', 'fourier', 'interpft' coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,.... ctrl_pfir_subband.nof_taps, ... - ctrl_pfir_subband.Wchan, ... + ctrl_pfir_subband.BWchan, ... ctrl_pfir_subband.r_pass, ... ctrl_pfir_subband.r_stop, ... ctrl_pfir_subband.coeff_w, ... @@ -77,93 +119,121 @@ coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,.... %ctrl_pfir_subband.coeff = flipud(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.gain = sum(ctrl_pfir_subband.coeff(:)) / ctrl_pfir_subband.downsample_factor; % Subband FFT parameters ctrl_pfft_subband.fft_size = tb.subband_fft_size; 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; -data_dt = []; -data_pfir_subband = []; -data_pfft_subband = []; -for block_nr = 1:tb.nof_blocks +data_dt = zeros(tb.nof_blocks, tb.subband_fft_size); +data_pfir_subband = zeros(tb.nof_blocks, tb.subband_fft_size); +data_pfft_subband = zeros(tb.nof_blocks, tb.nof_subbands); +for bi = 1:tb.nof_blocks % 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 end % Data path (DP) - [ctrl_wg, block_wg] = wg( ctrl_wg); - [ctrl_dt, block_dt] = dt( ctrl_dt, block_wg); + [ctrl_wg, block_wg] = wg( ctrl_wg); + [ctrl_dt, block_dt] = dt( ctrl_dt, block_wg); [ctrl_bsn_source] = bsn_source(ctrl_bsn_source); - [ctrl_pfir_subband, block_pfir_subband] = pfir(ctrl_pfir_subband, block_dt); - [ block_pfft_subband] = pfft(ctrl_pfft_subband, block_pfir_subband); + [ctrl_pfir_subband, block_pfir_subband] = pfir( ctrl_pfir_subband, block_dt); + [ block_pfft_subband] = pfft( ctrl_pfft_subband, block_pfir_subband); % Capture data at each DP interface - data_dt = [data_dt block_dt]; - data_pfir_subband = [data_pfir_subband block_pfir_subband]; - data_pfft_subband = [data_pfft_subband block_pfft_subband]; + data_dt(bi, :) = block_dt; + data_pfir_subband(bi, :) = block_pfir_subband; + data_pfft_subband(bi, :) = block_pfft_subband; end t_stop = cputime; disp(sprintf('Total processing time: %f seconds', t_stop-t_start)); % Plot data path results -ts = (0:tb.block_size*tb.nof_blocks-1)/tb.block_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 -tblock = (0:tb.nof_blocks-1); % time in blocks +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 +tblock = (0:tb.nof_blocks-1); % time in blocks % Plot DT output fig=fig+1; figure('position', [300+fig*20 200-fig*20 1000 800]); figure(fig); -plot(ts, data_dt) -ylim(1.2*[-ctrl_wg.ampl ctrl_wg.ampl]); -title('Delay tracking output data'); +data = data_dt'; +plot(ts, data(:)) +ylim([-1.3 1.3]); +title(sprintf('Delay tracking output data (WG subband %6.3f)', tb.subband_nr)); xlabel('Time [T bsn]'); ylabel('Voltage'); 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(ts, data_pfir_subband) -title('Subband polyphase FIR filter output - FFT input data'); +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; +figure('position', [300+fig*20 200-fig*20 1000 800]); +figure(fig); +data = data_pfir_subband'; +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]'); ylabel('Voltage'); 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; figure('position', [300+fig*20 200-fig*20 1000 800]); figure(fig); subplot(2,1,1); -plot(tsub, abs(data_pfft_subband)) +data = ampl'; +plot(tsub, data(:)) title('Subband data - amplitude'); xlabel('Time [T bsn]'); ylabel('Voltage'); grid on; subplot(2,1,2); -% force phase of too small signals to 0 -ampl = abs(data_pfft_subband); -noise = ampl<0.1*max(ampl); -phase = angle(data_pfft_subband)*180/pi; -phase(noise) = 0; -plot(tsub, phase) +data = phase'; +plot(tsub, data(:)) +ylim([-180 180]) title('Subband data - phase'); xlabel('Time [T bsn]'); ylabel('Phase [degrees]'); grid on; % 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; figure('position', [300+fig*20 200-fig*20 1000 800]); figure(fig); -nr = round(tb.subband_nr); -phase = angle(downsample(data_pfft_subband, tb.nof_subbands, nr))*180/pi; 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]'); ylabel('Phase [degrees]'); grid on; @@ -173,18 +243,17 @@ fig=fig+1; figure('position', [300+fig*20 200-fig*20 1000 800]); figure(fig); data = abs(data_pfft_subband); -data = data / max(data); -data = db(data); -x = floor(min(data(data~=-Inf))); -%x = -120; -%data(data<x) = x; -max(data) -min(data) -data = reshape(data, tb.nof_subbands, tb.nof_blocks); -color_range = size(colormap, 1); +data_abs_max = max(data(:)); +data = db(data); % no need to scale data, range is already normalized +x = ctrl_spec_subband.db_low; +%x = floor(min(data(data~=-Inf))) +data(data<x) = x; mymap = jet(-x); -%mymap(1,:) = [0 0 0]; +%mymap(1,:) = [0 0 0]; % force black for dB(0) colormap(mymap); -imagesc(data); +imagesc(data',[x 0]); colorbar; +title(sprintf('Subband spectogram (max value = %f = %f dB)', data_abs_max, db(data_abs_max))); +xlabel('Time [T bsn]'); +ylabel('Subband');