diff --git a/applications/lofar2/model/FIR_LPF_ApertifCF.m b/applications/lofar2/model/FIR_LPF_ApertifCF.m
new file mode 100644
index 0000000000000000000000000000000000000000..496056f22bab6621db33d61b9c3fe07a7336fa08
--- /dev/null
+++ b/applications/lofar2/model/FIR_LPF_ApertifCF.m
@@ -0,0 +1,124 @@
+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2016
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%-----------------------------------------------------------------------------
+% Author: R. de Wild, 2015 (Original)
+% E. Kooistra, 2016
+%
+% Purpose : Calculate FIR filter coefficients for Apertif channel
+% filterbank using fircls1.
+% Description :
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% workspace initiation
+clear all
+close all
+fig=0;
+
+% preset FIR-filter parameters
+K = 64; % Number of channels
+Taps = 8; % Number of taps, even 4 could be sufficient (erko)
+M = Taps*K; % Number of filter coefficients (taps)
+h_ix = 1:M; % coefficient index
+nbit = 9; % bits/coefficient = 1 sign + (nbit-1) mantissa
+fs = 781250; % sampling clock = 800.0 MHz / 1024
+f_nyq = fs/2; % Nyquist frequency
+f_pb = fs/K/2; % end of passband
+f_chan = fs/K;
+
+% normalized frequency axis & normalized amplitude axis
+w_pb = f_pb/f_nyq; % = f_chan
+
+% The fircls1 pass band deviation and stop band deviation are positive
+% fractions of 1. The ripple in dB needs to be divided by 10 to get from
+% dB = 10*log10() the log10() and divided by 2 because it is a relative
+% voltage value (instead of a relative power value).
+% When the pass band ripple in dB is positive then the linear factor is
+% sligthly larger than 1, so subtract 1 to get the pass band deviation.
+% The pass band deviation defines the +- fraction around 1. The stop band
+% 'ripple' in dB is negative and defines the attenuation. The stop band
+% deviation is the fraction of 1 that remains in the stop band. Therefore
+% both the pass band ripple and the stop band attenuation result both in
+% linear deviations that are close to 0.
+% The fircls1 description defines DEVP is the maximum passband deviation
+% or ripple (in linear units) from 1, and DEVS is the maximum stopband
+% deviation or ripple (in linear units) from 0.
+r_pb = 10^(0.025)-1; % passband ripple = 0.5 dB
+r_sb = 10^(-3.0); % stopband ripple = -60.0 dB
+
+% compute M filter-coefficients
+h = fircls1(M-1, w_pb, r_pb, r_sb);
+hfs_abs = abs(fftshift(fft(h ,M)));
+
+% compute quantized filter coefficients
+hmax = max(h);
+hq = double( uencode(h, nbit, hmax, 'signed') );
+hqfs_abs = abs(fftshift (fft( hq ,M))) * sum(h)/sum(hq); % normalize to response for DC
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Print FIR-filter parameters
+sprintf([ ...
+ ' fs = sample rate = %6.2f [Hz] \n' , ...
+ ' f_chan = channel bandwidth = %6.2f [Hz]\n' , ...
+ ' r_pb = pass band deviation = %9.7f = %6.2f dB_W ripple\n' , ...
+ ' r_sb = stop band deviation = %9.7f = %6.2f dB_W attenuation\n' , ...
+ ' M = number of coefficients = %d'], ...
+ fs, f_chan, r_pb, 20*log10(1+r_pb), r_sb, 20*log10(r_sb), M)
+
+% Save FIR-filter-coefficients
+% N.B. - read '.txt'-file with 'wordpad' or use \r\n !!!
+file_name = ['FIR_LPF_ApertifCF.txt'];
+FIRtable = [ h(h_ix) ] ;
+fid = fopen(file_name , 'w');
+fprintf(fid ,'%14.10f \n', FIRtable);
+fclose(fid);
+
+% Plot FIR-filter coefficients
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(hq);
+title(['FIR filter coefficients (', num2str(nbit), ' bits)'] );
+grid on;
+
+% Plot FIR-filter amplitude characteristic
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+f = (h_ix - M/2-1) * fs/M / f_chan;
+plot ( f, 20.0 * log10 ( hfs_abs ), 'r-', ...
+ f, 20.0 * log10 ( hqfs_abs ), 'k-') ;
+xlo = f(1);
+xhi = f(M);
+xlo = -3;
+xhi = 3;
+xlim( [xlo xhi] );
+grid on;
+zoom on;
+title('FIR filter transfer function');
+text( xlo, -5, ['sampling rate = ', num2str(fs), ' [Hz]']);
+text( xlo, -10, ['channel bandwidth = ', num2str(f_chan), ' [Hz]']);
+text( xlo, -15, ['number of taps = ', num2str(Taps)]);
+text( xlo, -20, ['number of channels = ', num2str(K)]);
+text( xlo, -25, ['number of coefficients = ', num2str(M)]);
+text( xlo, -30, ['number of bits/coefficient = ', num2str(nbit)]);
+xlabel('Frequency [channels]');
+ylabel('Power [dB]');
+print -dpng 'arts_fir';
diff --git a/applications/lofar2/model/pfir_coeff.m b/applications/lofar2/model/pfir_coeff.m
new file mode 100644
index 0000000000000000000000000000000000000000..aaad53b2b40ee985d793f23dff60e5a9a89676b3
--- /dev/null
+++ b/applications/lofar2/model/pfir_coeff.m
@@ -0,0 +1,145 @@
+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2016
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%-----------------------------------------------------------------------------
+% Author: W. Lubberhuizen, 2005
+% E. Kooistra, 2016
+%
+% pfir_coeff - compute polyphase FIR filter coefficients
+% Computes the coefficients of the prototype FIR filter for a polyphase
+% filterbank with N fft channels, and subfilter length L and channel
+% bandwidth BWchan. Typically BWchan = 1/N or somewhat more have -3dB
+% reponse between channels.
+%
+% The FIR coefficients can be caluculated with Matlab 'fir1' or with
+% Matlab 'fircls1' dependend on config.design.
+%
+% For Matlab fircls1() r_pass and r_stop are the maximum passband and
+% stop band deviations, respectively.
+% For Matlab fircls1() and Q>1, the filter coefficients are computed by
+% interpolating a shorter filter with an upsampling factor Q. This
+% reduces computation time substantially. N, L and Q must be powers of 2.
+%
+% The FIR coefficients are normalized to fit in range +-1. Given that
+% the impulse response has a sync pulse shape this implies that the
+% largest coefficient is h_fir_max = 1. Most energy is passed on via
+% the center tap and therefore the filter DC gain per poly phase
+% h_fir_dc_gain is close to 1. The quantization maps h_fir_max to
+% optional input argument coeff_q_max which defaults to 2^(nof_bits-1)-1.
+%
+% For nof_bits>1, the FIR filter coefficients are quantized with
+% nof_bits bits.
+%
+% The returned value is a coefficients vector of length NL, which also
+% represents the impulse response of prototype FIR filter.
+%
+% The run_pfir_coeff.m shows how this pfir_coeff.m can be used.
+function coeff = pfir_coeff(N, L, BWchan, r_pass, r_stop, nof_bits, config, coeff_q_max)
+
+% Default options
+if ~exist('coeff_q_max', 'var'); coeff_q_max = 2^(nof_bits-1)-1; end; % the quantization maps h_fir_max to coeff_q_max
+
+% requested filter length
+M2=N*L;
+if strcmp(config.design, 'fir1')
+ Q = 1;
+elseif M2<=1024
+ Q = 1;
+else
+ Q = L; % use interpolation to speed up calculation
+ %Q = 1; % no interpolation
+end
+% initial filter length
+M1=N*L/Q;
+
+% pass bandwidth
+w_pb = Q * BWchan;
+
+% compute initial filter
+if strcmp(config.design, 'fir1')
+ h_comp = fir1(M1-1, w_pb);
+ disp(sprintf('pfir_coeff : FIR coefficients calculated using Matlab fir1.\n'))
+else
+ if strcmp(config.design_flag, '')
+ h_comp = fircls1(M1-1, w_pb, r_pass, r_stop);
+ else
+ h_comp = fircls1(M1-1, w_pb, r_pass, r_stop, config.design_flag);
+ end
+ disp(sprintf(['pfir_coeff : FIR coefficients calculated using Matlab fircls1.\n', ...
+ 'pfir_coeff : . pass band deviation = %e\n', ...
+ 'pfir_coeff : . stop band deviation = %e\n'], ...
+ r_pass, r_stop))
+end
+
+if Q<=1
+ h_fir = h_comp;
+ disp(sprintf('pfir_coeff : FIR coefficients calculated directly with no interpolation.\n'))
+else
+ % 1a) Use Matlab fourier interpolation method
+ h_interpft = interpft(h_comp,length(h_comp)*Q);
+
+ % 1b) Use DIY fourier interpolation method
+ f1=fft(h_comp);
+ f2=zeros(1, M2);
+ % copy the lower frequency half.
+ n=0:M1/2;
+ f2(1+n)=f1(1+n);
+ % to make the impulse response symmetric in time,
+ % we need to apply a small phase correction,
+ % corresponding to a shift in the time domain.
+ f2(1+n)=f2(1+n).*exp(-sqrt(-1)*(Q-1)*pi*n/M2);
+ % recreate the upper part of the spectrum from the lower part
+ f2(M2-n)=conj(f2(2+n));
+ % back to time domain
+ h_fourier = real(ifft(f2));
+
+ % 2) Use resample interpolation method
+ h_resample = resample(h_comp, Q, 1);
+
+ % select FIR coefficients from either interpolation method (all are good, but resample causes peek gratings in stopband when N>64)
+ if strcmp(config.interpolate, 'resample')
+ h_fir = h_resample;
+ disp(sprintf('pfir_coeff : FIR coefficients calculated with resample interpolation.\n'));
+ elseif strcmp(config.interpolate, 'fourier')
+ h_fir = h_fourier;
+ disp(sprintf('pfir_coeff : FIR coefficients calculated with DIY Fourier interpolation.\n'));
+ else
+ h_fir = h_interpft;
+ disp(sprintf('pfir_coeff : FIR coefficients calculated with Matlab Fourier interpolation.\n'));
+ end
+end;
+
+% normalize the coefficients to fit in range +-1
+disp(sprintf('pfir_coeff : Normalizing the PFIR coefficients to fit in range +-1.'));
+h_fir = h_fir / max(h_fir);
+h_fir_max = max(h_fir); % note: is 1 due to the normalization
+h_fir_min = min(h_fir);
+h_fir_dc_gain = sum(h_fir) / N;
+disp(sprintf(' . Maximum coefficient value : %f', h_fir_max));
+disp(sprintf(' . Minimum coefficient value : %f', h_fir_min));
+disp(sprintf(' . DC gain per polyphase : %f\n', h_fir_dc_gain));
+
+% quantize the coeffients
+if nof_bits>0
+ disp(sprintf('pfir_coeff : Quantizing FIR coefficients nof_bits = %d\n', nof_bits));
+ h_fir = quantize(h_fir, h_fir_max, nof_bits, 'half_away', 'clip', coeff_q_max); % note: h_fir_max is 1 due to the normalization
+end;
+
+% output FIR coefficient
+coeff = h_fir;
\ No newline at end of file
diff --git a/applications/lofar2/model/pfir_coeff_adjust.m b/applications/lofar2/model/pfir_coeff_adjust.m
new file mode 100644
index 0000000000000000000000000000000000000000..99c0cbe0ef58cbd58daf96a42410d7daf77d2d00
--- /dev/null
+++ b/applications/lofar2/model/pfir_coeff_adjust.m
@@ -0,0 +1,170 @@
+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2016
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%-----------------------------------------------------------------------------
+% Author: E. Kooistra, 2016
+%
+% pfir_coeff_adjust - Compute and adjust polyphase FIR filter coefficients
+% Computes the coefficients using pfir_coeff.m and then adjusts the
+% coefficients dependent on:
+%
+% config.hp_factor : = 1, is initial or default value for BWchan
+% With config.hp_factor the BWchan can be set relative to fchan = 1/N.
+% Default BWchan = config.hp_factor/N = 1/N = fchan.
+%
+% config.hp_adjust :
+% when true then enable automatic channel half power bandwidth
+% adjustment of config.hp_factor by iteration until the
+% config.hp_factor has reached an accuracy of hp_resolution.
+% For floating point nof_bits=0 the config.hp_factor can achieve an
+% arbitrary accurate solution. For fixed point nof_bits>0 the solution
+% can only be as good as possible, so then decreasing hp_resolution
+% does not change the coeff anymore. The input BWchan is then ignored
+% and the initial BWchan is then set to config.hp_factor/N.
+%
+% config.dc_adjust :
+% when true then after quantization the FIR coefficients are adjusted
+% such that each of the N polyphase FIR sections with L taps each has
+% the same DC response. This then ensures that for DC input the output
+% of all N polyphases will also be DC. The FFT in a polyphase
+% filterbank will then transform this DC to bin 0, and there will be
+% no affect on the other bins.
+% The DC adjustment algorithm first determines the median of the sums
+% per polyphase. For each polyphase the L taps are then DC adjusted
+% to this dc_polyphases_median. If the DC adjustment is > L then first
+% an equal coarse adjustment is made to all L taps. The remining DC
+% adjustment is then < L and applied by adjusting the outer taps by 1.
+%
+% The DC adjustment must not cause coefficient overflow. If it does
+% then the script will rerun with coeff_q_max sufficiently smaller
+% than the default 2^(nof_bits-1)-1 until it succeeds.
+%
+% The run_pfir_coeff.m shows how this pfir_coeff_adjust.m and
+% pfir_coeff.m can be used.
+%
+% For more information see help of pfir_coeff.m.
+function coeff = pfir_coeff_adjust(N, L, BWchan, r_pass, r_stop, nof_bits, config, coeff_q_max)
+
+% Default options
+if ~exist('coeff_q_max', 'var'); coeff_q_max = 2^(nof_bits-1)-1; end; % the quantization maps h_fir_max to coeff_q_max
+
+while true
+ if config.hp_adjust==false
+ % Calculate the FIR coefficients using the specified BWchan
+ coeff = pfir_coeff(N, L, BWchan, r_pass, r_stop, nof_bits, config, coeff_q_max);
+
+ else
+ % Automatically adjust channel half power bandwidth, so need db(sqrt(2)) = -3 dB gain instead of half gain
+ hp_radix = 2; % choose 2 for binary search, or 10 for 'decimal' search
+ hp_step = 1/hp_radix;
+ hp_resolution = 0.000001;
+ while true
+ BWchan = config.hp_factor/N;
+ coeff = pfir_coeff(N, L, BWchan, r_pass, r_stop, nof_bits, config, coeff_q_max);
+ hf_abs = abs(fftshift(fft(coeff / sum(coeff), N*L)));
+ fi_0 = N*L/2+1; % frequency index of center of channel, so 0 Hz
+ fi_p = fi_0 + L/2; % frequency index of edge of channel at +f_chan/2
+ hf_abs_p = hf_abs(fi_p); % gain edge of channel at +f_chan/2
+ disp(sprintf('hp_factor = %10.8f, hf_abs_p = %10.8f, hp_step = %10.8f\n', config.hp_factor, hf_abs_p, hp_step));
+ if hf_abs_p < sqrt(0.5)
+ config.hp_factor = config.hp_factor + hp_step;
+ else
+ if hp_step <= hp_resolution, break; end
+ config.hp_factor = config.hp_factor - hp_step;
+ hp_step = hp_step/hp_radix;
+ config.hp_factor = config.hp_factor + hp_step;
+ end
+ end
+ end
+
+ if config.dc_adjust==true && nof_bits>0
+ % Automatically equalize the DC response of the polyphases to have flat DC level at FFT input
+ % Reshape the coefficients in N polyphases with L taps each.
+ % Work with q_full_scale quantized coefficients to have q_bit unit quantization.
+ q_bit = 1;
+ q_full_scale = 2^(nof_bits-1); % maximum amplitude, magnitude
+ q_coeff = round(q_full_scale * reshape(coeff, N, L).');
+ % The DC response is the sum of the FIR coefficients per polyphase.
+ dc_polyphases = sum(q_coeff);
+ % The median will be used as aim for the dc adjustment, because most polyphases are already
+ % close to the median. For larger N the first and last polyphases deviate the most.
+ dc_polyphases_median = median(dc_polyphases);
+ % DC adjust each polyphase to the DC median
+ for I = 1:N
+ polyphase = q_coeff(:,I);
+ dc_polyphase = sum(polyphase);
+ dc_adjust = dc_polyphase - dc_polyphases_median;
+ % First equally adjust all L coefficients by an equal amount of
+ % dc_offset >= 0
+ dc_offset = floor(abs(dc_adjust)/L);
+ % Then adjust the dc_remainder < L, over dc_adjust number of
+ % coefficients from the outer to the inner taps of the polyphase
+ % by 1
+ dc_remainder = abs(dc_adjust) - L*dc_offset;
+ % Treat dc_adjust < 0, > 0 separately, nothing to do when dc_adjust = 0
+ if dc_adjust > 0
+ polyphase = polyphase - dc_offset;
+ if dc_remainder > 0
+ for J = 1:dc_remainder
+ if mod(J,2)==0
+ polyphase(J) = polyphase(J) - 1;
+ else
+ polyphase(L+1-J) = polyphase(L+1-J) - 1;
+ end
+ end
+ end
+ end
+ if dc_adjust < 0
+ polyphase = polyphase + dc_offset;
+ if dc_remainder > 0
+ for J = 1:dc_remainder
+ if mod(J,2)==0
+ polyphase(J) = polyphase(J) + 1;
+ else
+ polyphase(L+1-J) = polyphase(L+1-J) + 1;
+ end
+ end
+ end
+ end
+ dc_polyphase = sum(polyphase);
+ dc_adjusted = dc_polyphase - dc_polyphases_median;
+ if dc_adjusted ~= 0
+ % Return with Error messages when dc_adjusted for a polyphase is not 0
+ disp(sprintf('Error polyphase %4d : dc_adjust = %4d, dc_offset= %4d, dc_remainder = %4d, dc_adjusted = %4d must be 0', I, dc_adjust, dc_offset, dc_remainder, dc_adjusted));
+ return
+ end
+ q_coeff(:,I) = polyphase / q_full_scale; % back to normalized full scale coefficients
+ end
+ coeff = reshape(q_coeff.', 1, N*L);
+
+ % Check whether the q_coeff are still in range q_min : q_max
+ q_full_scale = 2^(nof_bits-1);
+ q_max = q_full_scale-1;
+ q_margin = q_max - max(q_coeff(:));
+ if q_margin < 0
+ coeff_q_max = coeff_q_max+q_margin; % Recalculate pfir_coeff with more positive margin to fit DC adjustment
+ else
+ break; % Finished DC adjustment of the coefficients
+ end
+ else
+ break; % No need to while loop
+ end
+end
+
+end
diff --git a/applications/lofar2/model/readme_lofar2_model.txt b/applications/lofar2/model/readme_lofar2_model.txt
index 3dfaed1b1ef6aa4b4737e4f57987c5aef5e7e28e..82d595afef2aefb62f9cf0030e46a64045882bd2 100644
--- a/applications/lofar2/model/readme_lofar2_model.txt
+++ b/applications/lofar2/model/readme_lofar2_model.txt
@@ -1,5 +1,13 @@
This directory contains scripts and documents for modelling DSP aspects of LOFAR2.0:
+Contents:
+
+1) Original LOFAR1 scripts from FilterTaskForce.zip [1]
+2) Original LOFAR1 documents
+3) Copied from APERTIF [3]
+4) LOFAR2.0
+
+
References:
[1] https://git.astron.nl/desp/hdl/-/blob/master/applications/lofar1/FilterTaskForce.zip
@@ -9,6 +17,8 @@ References:
[4] https://git.astron.nl/desp/hdl/-/blob/master/libraries/dsp/verify_pfb/tb_tb_verify_pfb_wg.vhd
+
+
1) Original LOFAR1 scripts from FilterTaskForce.zip [1]
- git/hdl/libraries/dsp/filter/src/hex/Coeffs16384Kaiser-quant.dat :
@@ -32,14 +42,23 @@ References:
- quanreq.pdf, 2002 [2] : Modeling the Effects of Re-quantization for
ALMA-FC: Cascading of Degradation Factors
+
3) Copied from APERTIF [3]
APERTIF uses the same PFB FIR coefficients as LOFAR1.
+ FIR_LPF_ApertifCF.m Script Calculate FIR filter coefficients for Apertif channel filterbank
+ using fircls1 (original by Ronald de Wild).
+ pfir_coeff.m Function Compute polyphase filterbank coefficients (same purpose as
+ FIR_LPF_ApertifCF.m but with more options)
+ pfir_coeff_adjust.m Function Compute polyphase filterbank coefficients with optionbal half power
+ bandwidth adjustment (calls pfir_coeff.m)
pfir_coeff_dc_adjust.m Function Apply DC adjustment to polyphase filterbank coefficients
run_pfir_coeff.m Script Run pfir_coeff_adjust.m and pfir_coeff.m to show filter response in
time and frequency domain, and optionally calls pfir_coeff_dc_adjust.m
+
+
4) LOFAR2.0
- tb_tb_verify_pfb_wg.vhd, 2021 [4] : VHDL simulation of LOFAR1 PFB2 and APERTIF WPFB for LOFAR2.0