From ac7261b12dff8021b9b842bedc6f381de450305b Mon Sep 17 00:00:00 2001 From: Eric Kooistra <kooistra@astron.nl> Date: Wed, 5 Jan 2022 11:37:34 +0100 Subject: [PATCH] Copied FIR coef scripts from Apertif. --- applications/lofar2/model/FIR_LPF_ApertifCF.m | 124 +++++++++++++ applications/lofar2/model/pfir_coeff.m | 145 +++++++++++++++ applications/lofar2/model/pfir_coeff_adjust.m | 170 ++++++++++++++++++ .../lofar2/model/readme_lofar2_model.txt | 19 ++ 4 files changed, 458 insertions(+) create mode 100644 applications/lofar2/model/FIR_LPF_ApertifCF.m create mode 100644 applications/lofar2/model/pfir_coeff.m create mode 100644 applications/lofar2/model/pfir_coeff_adjust.m diff --git a/applications/lofar2/model/FIR_LPF_ApertifCF.m b/applications/lofar2/model/FIR_LPF_ApertifCF.m new file mode 100644 index 0000000000..496056f22b --- /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 0000000000..aaad53b2b4 --- /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 0000000000..99c0cbe0ef --- /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 3dfaed1b1e..82d595afef 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 -- GitLab