Use reorder.m to select the tb.subband_i.

applications/apertif/matlab/fringe_stopping.m

@@ -47,20 +47,22 @@ clear all;
 close all;
 fig=0;
 
+tb.model             = 'floating point';
+tb.model             = 'fixed point';
+
 tb.nof_complex       = 2;
 tb.nof_subbands      = 512;
 tb.bsn_init          = 0;
 tb.subband_nr        = 113;                             % subband range 0:511, can be fraction
-%tb.subband_nr        = 113.3;
+%tb.subband_nr        = 113.5;
+tb.subband_i         = round(tb.subband_nr);
 tb.subband_fft_size  = tb.nof_complex*tb.nof_subbands;  % subband filterbank real FFT
 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
-ctrl_wg.agwn_sigma   = 0;      % AGWN sigma
 if strcmp(tb.model, 'floating point')
     ctrl_wg.data_w            = 0;
+    ctrl_wg.agwn_sigma        = 0;      % AGWN sigma
     %ctrl_wg.agwn_sigma        = 0.01;
     ctrl_pfir_subband.coeff_w = 0;
     ctrl_pfir_subband.scale_w = 0;
@@ -71,8 +73,9 @@ if strcmp(tb.model, 'floating point')
 else
     ctrl_wg.data_w            = 8;
     lsb = 1/2^ctrl_wg.data_w;
-    %ctrl_wg.agwn_sigma        = 2.6*lsb;  % AGWN sigma
-    ctrl_pfir_subband.coeff_w = 9;      % bits/coefficient = 1 sign + (w-1) mantissa
+    ctrl_wg.agwn_sigma        = 0;      % AGWN sigma
+    %ctrl_wg.agwn_sigma        = 2.6*lsb;
+    ctrl_pfir_subband.coeff_w = 16;     % 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;
@@ -111,11 +114,22 @@ 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.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      = 'lofar file';
+ctrl_pfir_subband.config.design      = 'fir1';
+ctrl_pfir_subband.config.design      = 'fircls1';   % 'fir1', 'fircls1', 'lofar file'
 ctrl_pfir_subband.config.design_flag = 'trace';     % 'trace'
 ctrl_pfir_subband.config.interpolate = 'interpft';  % 'resample', 'fourier', 'interpft'
 
+if strcmp(ctrl_pfir_subband.config.design, 'lofar file')
+    % Load the quantized FIR coefficients of LOFAR station
+    ctrl_pfir_subband.downsample_factor  = 1024;
+    ctrl_pfir_subband.nof_taps           = 16;
+    ctrl_pfir_subband.nof_coefficients   = ctrl_pfir_subband.downsample_factor*ctrl_pfir_subband.nof_taps;
+    ctrl_pfir_subband.coeff_w            = 16;
+    coeff = load('data/Coeffs16384Kaiser-quant.dat');
+    coeff = reshape(coeff, ctrl_pfir_subband.downsample_factor, ctrl_pfir_subband.nof_taps);
+else
+    % Calculate FIR coefficients
     coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,....
                        ctrl_pfir_subband.nof_taps, ...
                        ctrl_pfir_subband.BWchan, ...
@@ -123,6 +137,7 @@ coeff = pfir_coeff(ctrl_pfir_subband.downsample_factor,....
                        ctrl_pfir_subband.r_stop, ...
                        ctrl_pfir_subband.coeff_w, ...
                        ctrl_pfir_subband.config);
+end
 
 %ctrl_pfir_subband.coeff = fliplr(coeff);
 %ctrl_pfir_subband.coeff = flipud(coeff);
@@ -135,11 +150,17 @@ ctrl_pfft_subband.fft_size = tb.subband_fft_size;
 ctrl_pfft_subband.complex  = false;
 ctrl_pfft_subband.gain = ctrl_pfft_subband.fft_size;
 
+% Reorder subband select parameters
+ctrl_reorder_subband.select = tb.subband_i + [1 2];
+ctrl_reorder_subband.select = tb.subband_i + 1;
+ctrl_reorder_subband.block_size = length(ctrl_reorder_subband.select);
+
 % Run the data path processing (one block per row)
 t_start = cputime;
 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);
+data_reorder_subband = zeros(tb.nof_blocks, ctrl_reorder_subband.block_size);
 for bi = 1:tb.nof_blocks
     % Control
     if ctrl_bsn_source.bsn == 30
@@ -152,11 +173,13 @@ for bi = 1:tb.nof_blocks
     [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_reorder_subband, block_reorder_subband] = reorder(   ctrl_reorder_subband, block_pfft_subband);
 
     % Capture data at each DP interface
     data_dt(bi, :)              = block_dt;
     data_pfir_subband(bi, :)    = block_pfir_subband;
     data_pfft_subband(bi, :)    = block_pfft_subband;
+    data_reorder_subband(bi, :) = block_reorder_subband;
 end
 t_stop = cputime;
 disp(sprintf('Total processing time: %f seconds', t_stop-t_start));
@@ -178,7 +201,7 @@ xlabel('Time [T bsn]');
 ylabel('Voltage');
 grid on;
 
-% Plot PFIR transfer function
+% Plot subband PFIR transfer function
 h = ctrl_pfir_subband.coeff;
 h = h(:);
 NL = ctrl_pfir_subband.nof_coefficients;
@@ -192,18 +215,18 @@ 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]);
+title(['Subband FIR filter transfer function for ', ctrl_pfir_subband.config.design]);
 xlabel('Frequency [channels]');
 ylabel('Magnitude [dB]');
 
-% Plot PFIR output
+% Plot subband 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]);
+ylim([-2 2]);             % DT step when tb.subband_nr is .5 causes double range
 xlabel('Time [T bsn]');
 ylabel('Voltage');
 grid on;
@@ -221,7 +244,7 @@ subplot(2,1,1);
 data = ampl';
 plot(tsub, data(:))
 title('Subband data - amplitude');
-xlabel('Time [T bsn]');
+xlabel(sprintf('Subband 0:%d at time [T bsn]', tb.nof_subbands-1));
 ylabel('Voltage');
 grid on;
 subplot(2,1,2);
@@ -229,20 +252,19 @@ data = phase';
 plot(tsub, data(:))
 ylim([-180 180])
 title('Subband data - phase');
-xlabel('Time [T bsn]');
+xlabel(sprintf('Subband 0:%d at time [T bsn]', tb.nof_subbands-1));
 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;
+phase = angle(data_pfft_subband(:, tb.subband_i + 1))*180/pi;
 
 fig=fig+1;
 figure('position', [300+fig*20 200-fig*20 1000 800]);
 figure(fig);
 plot(tblock, phase, '-o')
 ylim([-180 180])
-title(sprintf('Subband phase for subband %d in range 0:%d', nr, tb.nof_subbands-1));
+title(sprintf('Subband phase for subband %d in range 0:%d', tb.subband_i, tb.nof_subbands-1));
 xlabel('Time [T bsn]');
 ylabel('Phase [degrees]');
 grid on;