diff --git a/applications/lofar2/model/pfs_coeff_final.m b/applications/lofar2/model/pfs_coeff_final.m
index 7cc4798a14518588eb688864ce57a6dbe8f54791..3c3bb8947306f96fc3b736bdbed837ae719ad32c 100644
--- a/applications/lofar2/model/pfs_coeff_final.m
+++ b/applications/lofar2/model/pfs_coeff_final.m
@@ -1,128 +1,410 @@
+%------------------------------------------------------------------------------
+%
+% Copyright 2021
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+% http://www.apache.org/licenses/LICENSE-2.0
+%
+% Unless required by applicable law or agreed to in writing, software
+% distributed under the License is distributed on an "AS IS" BASIS,
+% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+% See the License for the specific language governing permissions and
+% limitations under the License.
+%
+%------------------------------------------------------------------------------
+%
+% Authors :
+% Created by FilterTaskForce team 2004 for LOFAR1 [1, 2]
+% Updated in 2022 for LOFAR20 by E. Kooistra
+% Purpose:
+%
+% References:
+% [1] hdl/applications/lofar1/FilterTaskForce/readme_FilterTaskForce.txt
+% [2] hdl/applications/lofar1/FilterTaskForce/pfs_coeff_final.m (Original
+% file from LOFAR1)
+% [3] LOFAR_MEM_162_digfilter_stopband3.pdf
+% [4] data/Coeffs16384Kaiser-quant.dat original LOFAR1 FIR coefficients
+% [5] data/Coefficient_16KKaiser.gold FTF reference FIR coefficients
+% [6] run_pfir_coeff.m, similar to this FIR filter generation script, used
+% for Apertif.
+%
+% Description:
+%
+% The LOFAR1 coefficients [4] were probably generated with this script, but
+% the exact same values cannot be reproduced, probably due to differences
+% in parameter settings. Nevertheless the FIR filter design method in this
+% scriptis still valid, so if we should need to generate new coefficients,
+% then we can do so with this script.
+%
+% Using reproduce_ftf = true does reproduce the FTF coefficients from [5].
+% However, these contain an error in r_pass, which could also be present the
+% LOFAR1 coefficients [4], which means that the r_pass ripple could be
+% improved.
+%
+% The half power subband bandwidth for the LOFAR1 [4] FIR filter is HP_BW =
+% 0.9954 * f_sub and for the FTF [5] FIR filter HP_BW = 0.9947 * f_sub.
+% Both FIR filters have the same maximum coefficient value of 2^15 - 1 =
+% 32767.
+%
+% With a plain FFT one expects a processing gain of N_point for the signal
+% to noise ration (SNR) of a bin, because the wanted signal is in that bin,
+% while the noise is still in all N_point bins. However, in poly phase
+% filterbank (PFB) simulations I (Eric) noticed that the PFB does not achieve
+% a processing gain of N_point in SNR. I think this is due to that the poly
+% phase FIR filter part in the PFB has a DC response that is not equal for
+% all N_point polyphases, that form the N_point input to the FFT part of the
+% PFB. There is a small variation in DC reponse that then causes leakage into
+% other bins and thus reduces the processing gain of the PFB to less than
+% N_points. This script shows that a more flat DC response over the poly
+% phases can be achieved by using a stronger window function for the FIR
+% coefficients calculation. However, more windowing would also increases the
+% half power bandwidth of the subbands. For LOFAR the reduced PFB processing
+% gain is no problem, because the PFB still meets the subband stop band
+% attenuation requirement of [3].
+%
+% The LOFAR1 FIR coefficients [4] have nof_bits = 16 bits. This is enough to
+% meet the specification of [3]. Therefore it is not necessary to regenerate
+% the FIR coefficients to 18 bit, even though the multipliers in the FPGA
+% could support 18b nearly for free.
+%
+% For comparision between LOFAR1 and LOFAR20 stations it is preferred to
+% keep using the same subband FIR coefficients. This is fine, because they
+% meet the specification of [3].
+%
+% Conclusion:
+% a) We can and should keep the original LOFAR1 subband FIR filter coefficients
+% also for LOFAR20.
+% b) We cannot exactly reproduce the LOFAR1 subband FIR filter coefficients, but
+% with this script we could generate an equivalent or improved set of subband
+% FIR coefficients if we would have to.
+%
+% Modifications:
+% a) Corrected missing -1 in r_pass=10^(r_pass_dB/20) - 1
+% b) Change r_stop_dB = -89 dB into -91 dB conform [3]
+% c) Clarified plots and corrected the f1_axis center
+% d) Support round() instead of uencode() for coefficient quantization
+% e) Compare created coefficients with LOFAR1 coefficients
+% f) Added plot of DC response per poly phase
+% g) Determine half power bandwidth of the pass band
+% h) Support loading and analysing LOFAR1 coefficients for comparison
+% i) Support changing the half power subband band width via relative_bw = 1
+
close all;
clear all;
+fig = 0;
+
+reproduce_ftf = true;
+%reproduce_ftf = false;
+
+%% User settings
+if not(reproduce_ftf)
+ % passband ripple (in dB);
+ r_pass_dB = 0.5;
+ r_pass = 10^(r_pass_dB/20) - 1;
+ % stopband ripple (in dB);
+ r_stop_dB = -89;
+ r_stop = 10^(r_stop_dB/20);
+ % relative subband band width
+ relative_bw = 1;
+ % word size for coefficients
+ nof_bits = 16;
+ % window function for coefficients
+ kaiser_beta = 1;
+ use_window = 'kaiser';
+ % quantization method
+ use_uencode = false;
+ % support analysing LOFAR1 coefficients
+ use_lofar1 = false;
+else
+ disp(sprintf('NOTE: recreate FTF coefficients result for reference'));
+ r_pass = 10^(0.5/20); % note missing -1
+ r_stop = 10^(-89/20); % -89 instead of -91 in [3]
+ relative_bw = 1;
+ nof_bits = 16;
+ kaiser_beta = 1;
+ use_window = 'kaiser';
+ use_uencode = true;
+ use_lofar1 = false;
+end
-fs =200*10^6;
+coeff_max = 2^(nof_bits - 1) - 1;
+
+
+%% Filter specification
% fft_size;
-N=1024;
+N = 1024;
% taps per frequency channel
-L=16;
-% passband ripple (in dB);
-r_pass=10^(0.5/20);
-% stopband ripple (in dB);
-r_stop=-89;
+L = 16;
% upsampling factor
-Q=16;
-% word size for coefficients
-W=16;
-
+Q = 16;
+% pass bandwidth
+subband_bw = relative_bw / N;
+comp_bw = Q * subband_bw;
% computed filter length
-M1=N*L/Q;
+M1 = N*L/Q;
% interpolated filter length
-M2=N*L;
-Nfft = M2;
-% compute filter
-c=fircls1(M1-1,Q/N,10^(r_pass/20),10^(r_stop/20),'trace');
-
-figure(30)
-Nfft =length(c);
-hold on
-spectrum_Fircls1= abs(fftshift(fft(c/max(c),Nfft)))/max(abs(fftshift(fft(c/max(c),Nfft))));
-plot([-Nfft/2+1:Nfft/2]/Nfft*fs,20*log10(spectrum_Fircls1),'r')
-hold off
-legend('Other candidate','filter given by fircls1')
-xlabel('f(Hz)')
-ylabel('Gain(dB)')
-title('Compirason of the prototype filter window before FFT interpolation');
-
-
-
-% use fourier interpolation to create final filter
-f1=fft(kaiser(length(c),1).'.*c);%hanning(length(c)).'.*
-f2=zeros(1,M2);
-% copy the lower frequency half.
-n=0:M1/2;
-
-figure(31)
-subplot(2,1,1)
-plot(-length(f1)/2+1:length(f1)/2,20*log10(fftshift(abs(f1))))
-xlabel('f(Hz)')
-title('Gain of the filter after FFT')
-ylabel('Gain(dB)')
-subplot(2,1,2)
-plot(-length(f1)/2+1:length(f1)/2,unwrap(angle(fftshift(f1))))
-title('Phase of the filter after FFT')
-xlabel('f(Hz)')
-ylabel('Gain(dB)')
-
-f2(1+n)=f1(1+n);
-% to make the impulse response exactly symmetric, we need to do a phase
-% correction
-% to make the impulse response exactly symmetric, we need to do a phase correction
-f2(1+n)=f2(1+n).*exp(-sqrt(-1)*(Q-1)*pi*n*1/M2);
+M2 = N*L;
+% fine frequency resolution
+P = 10;
+M3 = M2 * P;
+% frequency scale
+fs1 = N/Q;
+fs2 = N;
+f1_axis = [-M1/2:M1/2-1]/M1 * fs1; % for compute filter
+f2_axis = [-M2/2:M2/2-1]/M2 * fs2; % for interpolated filter
+f3_axis = [-M3/2:M3/2-1]/M3 * fs2; % for fine resolution spectrum
+
+%% Plot possible window functions
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot([kaiser(M1, 0.5), kaiser(M1, 1), kaiser(M1, 2), kaiser(M1, 4), hanning(M1), blackman(M1)])
+legend('kaiser beta = 0.5', 'kaiser beta = 1', 'kaiser beta = 2', 'kaiser beta = 4', 'hanning', 'blackman')
+title('Window functions');
+xlabel('Sample')
+ylabel('Gain')
+grid on;
+
+
+%% Compute filter
+h_comp = fircls1(M1-1, comp_bw, r_pass, r_stop, 'trace');
+
+fft_hcomp = fft(h_comp);
+fftshift_comp = fftshift(fft_hcomp);
+spectrum_comp = abs(fftshift_comp) / max(abs(fftshift_comp));
+spectrum_comp_dB = 20*log10(spectrum_comp);
+phase_comp = unwrap(angle(fftshift_comp));
+
+%% Compute windowed filter
+if strcmp(use_window, 'blackman')
+ disp(sprintf('NOTE: use blackman window'));
+ h_window = blackman(M1);
+elseif strcmp(use_window, 'hanning')
+ disp(sprintf('NOTE: use hanning window'));
+ h_window = hanning(M1);
+elseif strcmp(use_window, 'kaiser')
+ disp(sprintf('NOTE: use kaiser window (beta = %4.2f)', kaiser_beta));
+ h_window = kaiser(M1, kaiser_beta);
+else % default 'none'
+ disp(sprintf('NOTE: use no window'));
+ h_window = ones(M1, 1); % no window is rectangular window
+end
+h_wcomp = h_window' .* h_comp;
+
+%% Plot frequency response of compute filter
+fft_wcomp = fft(h_wcomp);
+fftshift_wcomp = fftshift(fft_wcomp);
+spectrum_wcomp = abs(fftshift_wcomp) / max(abs(fftshift_wcomp));
+spectrum_wcomp_dB = 20*log10(spectrum_wcomp);
+phase_wcomp = unwrap(angle(fftshift_wcomp));
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(f1_axis, spectrum_comp_dB, 'r', f1_axis, spectrum_wcomp_dB, 'g')
+legend('computed','windowed')
+title('Power spectrum of computed and windowed filter');
+xlabel('f [Hz]')
+ylabel('Gain [dB]')
+grid on;
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(f1_axis, phase_comp, 'r', f1_axis, phase_wcomp, 'g')
+legend('computed','windowed')
+title('Phase spectrum of computed and windowed filter'); % no difference
+xlabel('f [Hz]')
+ylabel('Phase [rad]')
+grid on;
+
+
+%% Use fourier interpolation to create final filter
+f2 = zeros(1, M2);
+% copy the lower frequency half of the computed filter into f2.
+n = 0:M1/2;
+f2(1+n) = fft_wcomp(1+n);
+
+% apply phase correction (= time shift) to make the impulse response exactly symmetric
+f2(1+n) = f2(1+n) .* exp(-sqrt(-1)*(Q-1)*pi*n*1/M2);
% create the upper part of the spectrum from the lower part
-f2(M2-n)=conj(f2(2+n));
-
-figure(32)
-subplot(2,1,1)
-plot(-length(f2)/2+1:length(f2)/2,20*log10(fftshift(abs(f2))))
-xlabel('f(Hz)')
-title('Gain of the filter after FFT')
-ylabel('Gain(dB)')
-subplot(2,1,2)
-hold on
-plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))))
-plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))),'r')
-hold off
-title('Phase of the filter after FFT')
-xlabel('f(Hz)')
-ylabel('Gain(dB)')
+f2(M2-n) = conj(f2(2+n));
% back to time domain
-c2=real(ifft(f2));
+h_fourier = real(ifft(f2));
+
+% use same scale
+h_fourier = h_fourier * coeff_max / max(h_fourier);
% quantize the coefficients
-c3=c2/max(c2);
-coeffs_quant = double(uencode(c3,16,max(c3),'signed'));
-% q=1/(2^(W-1));
-% c3=q*round(c3/q);
-% c3=c3*max(c2)/max(c3);
+if use_uencode
+ % uencode is equivalent to floor()
+ disp(sprintf('NOTE: use uencode() quantization'));
+ h_quant = double(uencode(h_fourier, nof_bits, coeff_max, 'signed'));
+else
+ % round is equivalent to floor(x + 0.5) provided that exactly 0.5 does not
+ % because round(0.5) = floor(x + 0.5) = 1, but round(-0.5) = -1.
+ disp(sprintf('NOTE: use round() quantization'));
+ h_quant = round(h_fourier);
+end
+
+%% Load LOFAR1 coefficients
+if use_lofar1
+ h_quant = load('data/Coeffs16384Kaiser-quant.dat');
+ h_quant = h_quant';
+end
+
+%% Check symmetry of FIR coefficients
+n=0:M2-1;
+h_reverse = h_quant(M2-n);
+d = h_quant == h_reverse;
+if sum(d) == length(d)
+ disp(sprintf('NOTE: quantized FIR coefficients are symmetrical.'));
+else
+ disp(sprintf('WARNING: quantized FIR coefficients are NOT symmetrical.'));
+end
+
+% save the coefficients
fid = fopen('data/Coefficient_16KKaiser.dat','w');
-fprintf(fid,'%i\n', coeffs_quant);
+fprintf(fid,'%i\n', h_quant);
fclose(fid);
-fprintf('\n\nWARNING !\n')
-fprintf('This coefficient recordered in Coefficient_16K.dat are scaled to one. It is necessary then to multiply them with the number of bits needed for the quanization: 2^N_bits \n\n')
+%% Compare with reference FTF coefficients
+h_ftf = load('data/Coefficient_16KKaiser.gold');
+d = h_ftf' == h_quant;
+if sum(d) == length(d)
+ disp(sprintf('NOTE: quantized FIR coefficients are identical to FTF coefficients.'));
+else
+ if reproduce_ftf
+ disp(sprintf('ERROR: quantized FIR coefficients are NOT identical to FTF coefficients.'));
+ else
+ disp(sprintf('NOTE: quantized FIR coefficients are NOT identical to FTF coefficients.'));
+ end
+end
+
+%% Compare with LOFAR1 coefficients
+h_lofar1 = load('data/Coeffs16384Kaiser-quant.dat');
+d = h_lofar1' - h_quant;
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(d)
+title('Difference between recreated coefficients and LOFAR1 coefficients');
+xlabel('Sample')
+ylabel('Coefficient')
+grid on;
-% freqz(c3,1,M2*32,200*10^6);
-% q=1/(2^(W-1));
-% c3=q*round(c3/q);
-% c3=c3*max(c2)/max(c3);
+d = h_lofar1' == h_quant;
+if sum(d) == length(d)
+ disp(sprintf('NOTE: quantized FIR coefficients are identical to LOFAR1 coefficients.'));
+else
+ disp(sprintf('NOTE: quantized FIR coefficients are NOT identical to LOFAR1 coefficients.'));
+end
-figure(2);
-n=0:M2-1;
-plot([c2',c2(M2-n)',c3']);
-title('impulse response');
-legend('unquantized','reversed','quantized');
+
+%% Plot (un)quantized impulse response
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot([h_fourier', h_quant']);
+legend('unquantized','quantized');
+title('Impulse response');
+xlabel('Sample')
+ylabel('Coefficient')
+grid on;
+
+%% Plot quantization error histogram to see effect of use_uencode
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+hist(h_fourier - h_quant,100);
+title('Quantization error histogram');
+xlabel('Quantization error');
+ylabel('Relative probability');
+grid on;
+
+
+%% Plot DC response per polyphase
+% . the DC response over the polyphases gets more flat with more windowing
+dc_fourier = sum(reshape(h_fourier, N, L)');
+dc_quant = sum(reshape(h_quant, N, L)');
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot([dc_fourier', dc_quant']);
+legend('unquantized','quantized');
+title('DC response per polyphase');
+xlabel('Polyphase index')
+ylabel('DC response')
grid on;
-figure(2);
-hist(c2-c3,100);
-title('quantization error histogram');
-xlabel('quantization error');
-ylabel('relative probability');
-
-figure(10);
-freqz(c3,1,M2*4);
-title('frequency response');
-Nfft = 16384*8;
-c3 = c3/max(c3);
-%spectrum_fft = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
-spectrum = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
-frequency = [-Nfft/2+1:Nfft/2]*200*10^6/Nfft;
-plot(frequency(length(frequency)/2:end),20*log10(spectrum(length(frequency)/2:end)))
-grid on
-title('Gain of the filter window with quantized coefficients (16bits) using Fircls1 window with the FFT interpolation method: blackman window')
-xlabel('f(Hz)')
-ylabel('Gain(dB)')
-legend('FFT interpolation')
+
+%% Plot frequency response of original and quantized filter
+
+%freqz(h_fourier, 1, M3);
+%freqz(h_quant, 1, M3);
+
+fft_org = fft(h_fourier, M3);
+fftshift_org = fftshift(fft_org);
+spectrum_org = abs(fftshift_org) / max(abs(fftshift_org));
+spectrum_org_dB = 20*log10(spectrum_org);
+phase_org = unwrap(angle(fftshift_org));
+
+fft_quant = fft(h_quant, M3);
+fftshift_quant = fftshift(fft_quant);
+spectrum_quant = abs(fftshift_quant) / max(abs(fftshift_quant));
+spectrum_quant_dB = 20*log10(spectrum_quant);
+phase_quant = unwrap(angle(fftshift_quant));
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(f3_axis, spectrum_org_dB, 'r', f3_axis, spectrum_quant_dB, 'g')
+legend('original','quantized')
+title('Power spectrum of original and quantized filter');
+xlabel('f [Hz]')
+ylabel('Gain [dB]')
+grid on;
+
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+% There is a large absolute phase offset between phase_org and phase_quant,
+% due to the quantization noise. Remove this absolute phase offset to be
+% able to compare the phase within the pass band, to see that it is the
+% same and is linear from -0.5:0.5. Use the phase offset at DC so the center
+% at index M3/2. Manually zoom in to range -1:1 to see the linear phase from
+% -0.5:0.5 in the subband.
+phase_offset = phase_org(M3/2) - phase_quant(M3/2);
+phase_quant = phase_quant + phase_offset;
+plot(f3_axis, phase_org, 'r', f3_axis, phase_quant, 'g')
+legend('original','quantized')
+title('Phase spectrum of original and quantized filter');
+xlabel('f [Hz]')
+ylabel('Phase [rad]')
+grid on;
+
+
+%% Determine half power bandwidth of the subband
+pass_band = spectrum_quant >= 0.5;
+pass_indices = find(pass_band);
+pass_index_lo = pass_indices(1);
+pass_index_hi = pass_indices(end);
+% linear interpolate to relative index at 0.5 power
+a = spectrum_quant(pass_index_lo);
+b = spectrum_quant(pass_index_lo - 1);
+pass_index_lo = pass_index_lo + (0.5 - b) / (a - b);
+a = spectrum_quant(pass_index_hi);
+b = spectrum_quant(pass_index_hi + 1);
+pass_index_hi = pass_index_hi + 1 - (0.5 - b) / (a - b);
+pass_bw = (pass_index_hi - pass_index_lo + 1) / M3 * fs2;
+%disp(sprintf('%f', pass_index_lo));
+%disp(sprintf('%f', pass_index_hi));
+disp(sprintf('NOTE: quantized FIR filter normalized subband bandwidth is %6.4f * f_sub', pass_bw));