Select Git revision
pfir_coeff_adjust.m
-
Eric Kooistra authoredEric Kooistra authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
pfir_coeff_adjust.m 8.22 KiB
%-----------------------------------------------------------------------------
%
% 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