Skip to content
Snippets Groups Projects
Commit ac7261b1 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Copied FIR coef scripts from Apertif.

parent c5602f5f
Branches
Tags
1 merge request!188Updated model/pfs_coeff_final.m, but preserved old FTF results for reference...
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/FIR_LPF_ApertifCF.m 0 → 100644
+ 124
− 0
View file @ ac7261b1
%-----------------------------------------------------------------------------
%
% 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';
applications/lofar2/model/pfir_coeff.m 0 → 100644
+ 145
− 0
View file @ ac7261b1
%-----------------------------------------------------------------------------
%
% 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
applications/lofar2/model/pfir_coeff_adjust.m 0 → 100644
+ 170
− 0
View file @ ac7261b1
%-----------------------------------------------------------------------------
%
% 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
applications/lofar2/model/readme_lofar2_model.txt
+ 19
− 0
View file @ ac7261b1
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment