Prepared working with Lofar1 pfs_coeff_final.m and Apertif run_pfir_coeff.m in lofar2/model dir.

applications/lofar1/FilterTaskForce/readme_FilterTaskForce.txt

+2021-12-7 Eric Kooistra (Jira L2SDP-215)
+
+Invertigation of FilterTaskForce.zip contents.
+
+1) Scope
+2) FilterTaskForce
+   a) Powerpoint presentations:
+   b) Word documents:
+   c) MATLAB files:
+   d) Coefficient files
+
+
+
+1) Scope
+
+I checked the FilterTaskForce.zip that Andre Gunst still had with documents
+and Matlab code for the Polyphase Filterbank that is used in LOFAR1.
+
+Unfortunately the documents and note are not very clear. The most useful
+things are:
+- the PFB drawings in filterbanks2.ppt
+- the filter specification in MeasurementsEffect.ppt
+
+The pfs_coeff_final.m generates the Coefficient_16KKaiser.dat, these FIR
+coefficients are almost the same as in:
+
+  git/hdl/libraries/dsp/filter/src/hex/Coeffs16384Kaiser-quant.dat
+
+that is used in LOFAR1, but unfortunately not exactly the same. Hence the
+most useful files are pfs_coeff_final.m and Coefficient_16KKaiser.dat. For
+further analysis in LOFAR2 these files are copied to:
+
+  git/hdl/applications/lofar2/model
+
+
+
+2a) FTF Powerpoint presentations:
+
+* Feb  9  2004 'Filter test bench algorithm verification.ppt' : not clear and not useful.
+* Feb 10  2004 '1KHz RFI Impact.ppt':
+
+  Plot of output power and phase over PFB bands, the phase is linear per band
+  and over all bands.
+
+* Feb 25  2004 'RFI study.ppt' : not clear and not useful.
+* Feb 27  2004  Filter_presentation.ppt : not clear and not useful.
+* May 14  2004 'Filter Task Force(final).ppt' : not clear and not useful.
+* May 14  2004  filterbanks2.ppt
+  Diagram of FIR and FFT for PFB, with L = M*N = 256 * 16 = 4096 coefficients.
+  . M polyphases and M point FFT
+  . N taps per polyhase
+
+* Jul  9  2004  FIRprototype.ppt : not clear and not useful.
+* Jul  9  2004  FIRprototype2.ppt : not clear and not useful.
+* Sep 17  2004  Filtermeeting2.ppt : not clear and not useful, except:
+
+  The filter design method fircls1 LSQ based methods gives good result for the
+  basic 2-4K design. An interpolation should be applied in order to realize
+  4K -> 16K coefficient generation. (Add 16K coefficients within the 4K
+  coefficients). Two methods are possible:
+  - Classic: 0 padding + FIR (long if good filter quality)
+  - FFT : creation of the bandwidth in padding zeros directly in the frequency domain.
+  ==> FFT interpolation yields best results
+
+* Oct 11  2004  MeasurementsEffect.ppt: not clear and not useful, except:
+
+  Stop band: -90.5dB
+  Pass band ripple: -0.45 dB
+  -3dB -> 40dB: 23.3%
+  -3dB -> 90.5dB: 32.89%
+
+* Feb 10  2005 'Quantization results 200MHz.ppt' : not useful.
+
+
+2b) FTF  Word documents:
+
+* Dec 15  2003  Memoonbasebandshift.doc
+  Idea to shift the bands per station, so that the subband aliasing decorrelates at
+  the centrals processing. Not used.
+
+* May 14  2004 'Problematic choice of the filter.doc'
+  Some knowns pros and cons of critically sampled and oversampled PFB
+
+* Nov 23  2004  MinuteFTF2806.doc : not clear and not useful.
+
+RTF documents : all not useful
+
+* Oct 27  2003 'minutes 3.rtf'
+* Jan 22  2004 'Minutes filter-22-01.rtf'
+* Feb 27  2004  minutesFTF26-02.rtf
+* Mar 10  2004 'RFI discussion.rtf'
+* Mar 22  2004 'FTF - minutes_22-03-2004.rtf'
+* Jun 28  2004 'Document Plan.rtf'
+
+Oct  7  2003 'FTFarchive/Lofar station digital filtering.doc' : Forming the Filter task force
+Nov 14  2003  FTFarchive/FIRresult14-11.doc : not useful.
+Dec  9  2003  FTFarchive/STW-09-12.doc : Tasks and persons overview
+
+Excel file :  all not useful
+
+* Oct 26  2004 filter2.xls
+* Oct 26  2004 filter.xls
+
+
+2c) FTF  MATLAB files:
+
+* Aug 16  2004 pfs_coeff.m : Blackman
+* Aug 19  2004 pfs_coeff2.m : Hanning
+* Oct 26  2004 pfs_coeff_final.m : Kaiser
+  - copied to: /dop466_0/kooistra/git/hdl/applications/lofar2/model/pfs_coeff_final.m
+  - generates Coefficient_16KKaiser.dat
+
+
+2d) FTF  Coefficient files
+
+* Aug 19  2004 Coefficient_16K.dat
+* Oct 27  2004 Coefficient_16KKaiser.dat
+  - copied to: /dop466_0/kooistra/git/hdl/applications/lofar2/model/data/Coefficient_16KKaiser.dat
+
+

applications/lofar2/model/pfir_coeff_dc_adjust.m

+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2018
+% 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, Apr 2018
+%
+% pfir_coeff_dc_adjust - Adjust the DC response per poly phase of an
+%   existing set of PFIR coefficients
+%
+%   The input coeffIn is an array of N*L integer PFIR coefficients that fit 
+%   in nof_bits. The number of polyphase is N and the number of taps per
+%   polyphase is L. The order of the input coeffIn is first all L
+%   coefficients for phase 0 then for phase 1, and last for phase N-1.
+%
+%   The coeffIn values 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 adjusted FIR coefficients are returned via coeffNoDC. The DC
+%   adjustment must not cause coefficient overflow. The nof_bits
+%   defines the maximum integer coefficient range. If the coeffNoDC does
+%   not fit, than the function reduces the input coeffIn by a factor
+%   and then recalculates coeffNoDC, until the coefficients all fit within
+%   nof_bits.
+%
+%   The run_pfir_coeff.m shows how this pfir_coeff_adjust.m can be used.
+function coeffNoDC = pfir_coeff_dc_adjust(coeffIn, N, L, nof_bits)
+
+iteration = 0;
+while true
+    iteration = iteration + 1;
+    % 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.
+    q_coeff = reshape(coeffIn, 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;
+    end
+    coeffNoDC = 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_no_dc_max = max(coeffNoDC);
+    if q_no_dc_max > q_max
+        coeffIn = round(coeffIn * q_max / q_no_dc_max);  % Scale coeff with more positive margin to fit DC adjustment
+    else
+        break;  % Finished DC adjustment of the coefficients
+    end
+end
+disp(sprintf('DC adjusted in %d iteration(s)', iteration));
+
+end

applications/lofar2/model/pfs_coeff_final.m

@@ -87,7 +87,7 @@ coeffs_quant = double(uencode(c3,16,max(c3),'signed'));
 % c3=q*round(c3/q);
 % c3=c3*max(c2)/max(c3);
 
-fid = fopen('Coefficient_16KKaiser.dat','w');
+fid = fopen('data/Coefficient_16KKaiser.dat','w');
 fprintf(fid,'%i\n', coeffs_quant);
 fclose(fid);
 fprintf('\n\nWARNING !\n')

applications/lofar2/model/readme_lofar2_model.txt

@@ -2,32 +2,44 @@ This directory contains scripts and documents for modelling DSP aspects of LOFAR
 
 References:
 
-[1] = found by Andre Gunst from LOFAR1
-[2] = https://svn.astron.nl/UniBoard_FP7/RadioHDL/trunk/applications/apertif/matlab/apertif_matlab_readme.txt
-[3] = https://support.astron.nl/jira/browse/L2SDP-212
-[4] = https://git.astron.nl/desp/hdl/-/blob/master/libraries/dsp/verify_pfb/tb_tb_verify_pfb_wg.vhd
+[1] https://git.astron.nl/desp/hdl/-/blob/master/applications/lofar1/FilterTaskForce.zip
+    found by Andre Gunst from LOFAR1,for overview see readme_FilterTaskForce.txt
+[2] from Andre Gunst
+[3] https://git.astron.nl/desp/hdl/-/blob/master/apertif_matlab/README.md
 
+[4] https://git.astron.nl/desp/hdl/-/blob/master/libraries/dsp/verify_pfb/tb_tb_verify_pfb_wg.vhd
 
-1) Original LOFAR1 scripts
+1) Original LOFAR1 scripts from FilterTaskForce.zip [1]
 
-- git/hdl/libraries/dsp/filter/src/hex/Coeffs16384Kaiser-quant.dat :  Original PFIR filter coefficients from LOFAR 1.0 subband filterbank
-                                                                      (16 tap, 1024 point FFT), the script that created these PFIR
-                                                                      coefficients is not available anymore (For Apertif Erko looked in old
-                                                                      repositories back to 2005 of Lofar software and Station firmware
-                                                                      but could not find it.)
+- git/hdl/libraries/dsp/filter/src/hex/Coeffs16384Kaiser-quant.dat :
+     Original PFIR filter coefficients from LOFAR 1.0 subband filterbank
+     (16 tap, 1024 point FFT), the script that created these PFIR
+     coefficients is not available anymore (For Apertif Erko looked in old
+     repositories back to 2005 of Lofar software and Station firmware
+     but could not find it.)
                                                                       
-- pfs_coeff_final.m   [1] : creates Coefficient_16KKaiser.dat, but using e.g. meld shows that these do differ slightly from
-                           Coeffs16384Kaiser-quant.dat that are actually used in LOFAR1.
+- pfs_coeff_final.m [1]
+     creates Coefficient_16KKaiser.dat, but using e.g. meld shows that
+     these do differ slightly within +-40 from
+     Coeffs16384Kaiser-quant.dat that are actually used in LOFAR1.
 
 
 
 2) Original LOFAR1 documents
 
-- lowbitfir.pdf, 2001 [1] : Simulation Results of Low Bit FIR Filters for ALMA
-- quanreq.pdf, 2002   [1] : Modeling the Effects of Re-quantization for ALMA-FC: Cascading of Degradation Factors
- 
- 
-3) LOFAR2.0
+- filterbanks2.ppt    [1] : With diagrams of PFB structure
+- lowbitfir.pdf, 2001 [2] : Simulation Results of Low Bit FIR Filters for ALMA
+- 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.
+
+  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
 
-- Jira story L2SDP-212 [3]          : Investigate SDP firmware quantization 
 - tb_tb_verify_pfb_wg.vhd, 2021 [4] : VHDL simulation of LOFAR1 PFB2 and APERTIF WPFB for LOFAR2.0