diff --git a/applications/apertif/matlab/run_pfb.m b/applications/apertif/matlab/run_pfb.m
index 2e063cdbe3542997a53e295522ac5f1a32d64f80..2761de36597a2e78868f9942d7faa295ea2ebb30 100644
--- a/applications/apertif/matlab/run_pfb.m
+++ b/applications/apertif/matlab/run_pfb.m
@@ -64,6 +64,7 @@ tb.subband_wg        = 1;           % subband range 0:tb.nof_subbands-1, can be
 tb.subband_wg        = 301;
 tb.subband_wg        = 301+1/13;
 tb.subband_wg        = 64;
+tb.subband_wg        = 64+1/64;
 
 tb.phase             = 0;           % initial 'sinusoid' phase offset normalized to 2pi
 tb.sop               = 1;           % initial 'impulse' start index in range ctrl_wg.block_size = tb.subband_fft_size
@@ -331,6 +332,40 @@ if strcmp(tb.model_quantization, 'fixed point')
     fclose(fid);
 end
 
+%% Log subband statistics for subband_i
+% bin indices for a subband spectrum after the impulse response has
+% settled, so after ctrl_pfir_subband.nof_taps * tb.nof_subbands
+tg = (ctrl_pfir_subband.nof_taps + 4) * tb.nof_subbands;
+tgI = tg + 1;
+sub_bins = tgI + [0:tb.nof_subbands-1];   % Matlab indices of subbands in block tg
+sub_I_bin = tgI + tb.subband_i;           % Matlab index of subband_I
+
+sub_ampl = abs(data_pfft_subband);
+data = sub_ampl.';
+data = data(:);
+sub_I_ampl = data(sub_I_bin);
+sub_I_power = sub_I_ampl^2;
+sub_spurious = data(sub_bins);
+sub_spurious(tb.subband_I) = 0;
+sub_spurious_max = max(sub_spurious);
+sub_spurious_power = sum(sub_spurious.^2);
+sub_I_snr = sub_I_power / sub_spurious_power;
+sub_I_snr_max = sub_I_power / sub_spurious_max^2;
+
+disp(sprintf(['\nStatistics for WG subband %f in subband bin %d:\n', ...
+              '. Subband amplitude            = %10.8f\n', ...
+              '. Subband power                = %10.8f\n', ...
+              '. Spurious amplitude max       = %10.8f\n', ...
+              '. Spurious total power         = %10.8f\n', ...
+              '. SNR subband / spurious total = %f dB\n', ...
+              '. SNR subband / spurious max   = %f dB\n'], ...
+                  tb.subband_wg, tb.subband_i, ...
+                  sub_I_ampl, ...
+                  sub_I_power, ...
+                  sub_spurious_max, ...
+                  sub_spurious_power, ...
+                  10 * log10(sub_I_snr), ...
+                  10 * log10(sub_I_snr_max)));
 
 %% Plot data path results
 xfig = 300;
@@ -346,6 +381,7 @@ tsub_one = (0:                    tb.nof_tsub-1);                       % time i
 sub_I      = tb.subband_I      + [0: tb.nof_subbands: tb.nof_subbands*tb.nof_tsub-1];  % get Matlab indices of all subband_I
 sub_Iplus1 = tb.subband_Iplus1 + [0: tb.nof_subbands: tb.nof_subbands*tb.nof_tsub-1];  % get Matlab indices of all next neighbour subband_I+1
 
+
 %% Plot WG output
 fig=fig+1;
 figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);