pfir_coeff_adjust.m

%-----------------------------------------------------------------------------
%
% 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