function b = firhalfband(N,fp,varargin)
%FIRHALFBAND  Halfband FIR filter design.

%   References:
%     [1] T. Saramaki, Finite impulse response filter design, in
%          Handbook for Digital Signal Processing. S.K. Mitra and
%          J.F. Kaiser Eds. Wiley-Interscience, N.Y., 1993, Chapter 4.

%   Copyright 1999-2023 The MathWorks, Inc.

narginchk(2,5);

flags = validateNparseinput(N,fp,varargin{:});

if flags.winflag
    % Design using the window method
    b = windowdesign(N,fp,flags,varargin{:}); % fp is window vector
elseif flags.firceqripflag
    % Design using firceqrip
    b = firceqripdesign(N,fp,flags); % fp is DEV
else
    % Design via firpm
    b = firpmdesign(N,fp,flags,varargin{:}); % fp is passband frequency
end

% Convert to highpass if so wanted
if ~flags.lowpassflag
    b = firlp2hp(b);
end

%-------------------------------------------------------------
function b = windowdesign(N,fp,flags,dev,varargin)
% Design via window method, don't use fir1 so we get exact zeros
% every other sample, and exactly 0.5 at the middle sample

if flags.minphase
    % Minphase option not supported for window method
    error(message('dsp:firhalfband:minphaseWin'));
elseif flags.zerophase
    % Nonnegative zero-phase option not supported for window method
    error(message('dsp:firhalfband:zerophaseWin'));
end

if flags.minordflag
    % Estimate the required order and design the filter using Kaiser window
    [N,Win] = estimatekaiserord(fp,dev);
    b = lowpasslband(N,2,Win);
    tw = 2*(0.5-fp);   
    Astop = -db(dev);  

    % Check if the designed filter meets the required Stopband Attenuation
    [~, ~, ~, ~, ~, Astop_m] = dsp.filterdesign.internal.lowpassmeasure(b, tw, 2,"Stopband");
    if Astop_m < Astop
        % Design the filter iteratively if condition on stopband
        % attenuation is not met
        b = dsp.filterdesign.internal.designNyquistMinOrder(2,tw,Astop);
        % Remove first and last coefficients if they are zeros
        if ~rem((length(b)-1)/2,2)
            b = b(2:end-1);
        end
    end
else
    b = lowpasslband(N,2,fp);
end


%-------------------------------------------------------------
function [N,fp] = estimatekaiserord(fp,Dev)
% Estimate minimum order filter to meet peak ripple


[N,~,beta] = kaiserord([fp,1-fp],[1 0],[Dev,Dev]);
% Increase N to next integer such that N/2 is odd integer
N = adjustord(N);

fp = kaiser(N+1,beta); % fp is a window vector in windowdesign

%-------------------------------------------------------------
function ddev = computeddev(dev)
ddev = dev^2/2*(0.5+dev^2/2)/0.5;

%-------------------------------------------------------------
function b = computeminphase(bnonneg,ddev)

N = ceil(length(bnonneg)/2);
bnonneg(N) = bnonneg(N) + ddev; % Make zerophase nonnegative
bnonneg = 0.5/(0.5+ddev)*bnonneg; % Normalize to have gain 0.5 at 0.5
b = firminphase(bnonneg);

%-------------------------------------------------------------
function b = firceqripdesign(N,dev,flags)
% Halfband design via firceqrip

ddev = computeddev(dev);
if flags.minphase
    bnonneg = firceqrip(2*N,0.5,[ddev,ddev]); % Design double-order filter
    b = computeminphase(bnonneg,ddev);
else
    b = firceqrip(N,0.5,[dev,dev]);
    % Force exact zeros and exact 0.5 in the middle
    b = forcehalfband(b);
end
%-------------------------------------------------------------
function b = firpmdesign(N,fp,flags,dev,varargin)
% Halfband design via firpm

if flags.minphase
    if flags.minordflag
        ddev = dev^2/2*(0.5+dev^2/2)/0.5;
        bnonneg = firgr('mineven',[0, fp, 1-fp, 1],[1 1 0 0],[ddev,ddev]);
        b = computeminphase(bnonneg,ddev);
    else
        b = minphasefirpmdesign(N,fp);
    end
elseif flags.zerophase
    b = zerophasefirpmdesign(N,fp);
else
    if flags.minordflag
        % Estimate the required order
        N = estimatefirpmord(fp,dev);
    end

    % If order is multiple of 4, design for order-2 and add 2 zeros
    mult4 = false;
    if rem(N,4) == 0
        mult4 = true;
        N = N - 2;
    end

    % Design a type 2 filter with one band
    btemplate = firgr(N/2,[0 2*fp],[.5 .5],'2');

    % Insert zeros every other sample
    b(1:2:N+1) = btemplate;

    % Add a discrete impulse of 1/2 height
    b(N/2+1) = 0.5;

    % Add two zeros for orders that are multiples of 4
    if mult4
        b = [0, b, 0];
    end
end
%-------------------------------------------------------------
function N = estimatefirpmord(fp,Dev)
% Estimate minimum order filter to meet peak ripple

N = firpmord([fp,1-fp],[1 0],[Dev,Dev]);

%Increase N to next integer such that N/2 is odd integer
N = adjustord(N);


%-------------------------------------------------------------
function b = minphasefirpmdesign(N,fp)

fedges = [fp, 1-fp];
[bnonneg,err] = firgr(2*N,[0 fedges 1],[1 1 0 0]); % Design double-order filter
b = computeminphase(bnonneg,err);

%-------------------------------------------------------------
function b = zerophasefirpmdesign(N,fp)

fedges = [fp, 1-fp];
[bnonneg,err] = firgr(N,[0 fedges 1],[1 1 0 0]); 
N = ceil(length(bnonneg)/2);
bnonneg(N) = bnonneg(N) + err; % Make zerophase nonnegative
b = 0.5/(0.5+err)*bnonneg; % Normalize to have gain 0.5 at 0.5
% Force exact zeros and exact 0.5 in the middle
b = forcehalfband(b);

%-------------------------------------------------------------
function b = forcehalfband(b)
% Force exact zeros and exact 0.5 in the middle
indx = ceil(length(b)/2);
b(indx+2:2:end) = 0;
b(indx-2:-2:1) = 0;
b(indx) = 0.5;

%-------------------------------------------------------------
function N = adjustord(N)
% Increase N to next integer such that N/2 is odd integer

if N < 6
    N = 6;
else
    while (N/2 ~= round(N/2)) || ~rem(N/2,2)
        N = N + 1;
    end
end
%-------------------------------------------------------------
function flags = validateNparseinput(N,fp,varargin)

% Generate defaults
flags.winflag = 0;
flags.minordflag = 0;
flags.kaiserflag = 0;
flags.lowpassflag = 1;
flags.firceqripflag = 0;
flags.minphase = 0;
flags.zerophase = 0;

% Check if highpass and/or minphase/nonnegative is wanted
if nargin > 2 && dsp.util.isChar(varargin{end})
    [flags,flagschanged] = checklaststr(flags,char(varargin{end}));
    if flagschanged
        varargin = varargin(1:end-1);
    end
end
% Check AGAIN if highpass and/or minphase/nonnegative is wanted
if 2 + length(varargin) > 2 && dsp.util.isChar(varargin{end})
    [flags,flagschanged] = checklaststr(flags,char(varargin{end}));
    if flagschanged
        varargin = varargin(1:end-1);
    end
end

% At this point the possible trailing 'low','high', 'minphase' or
% 'nonnegative' flags have been removed
switch 2+length(varargin)
    case 2
        % (N,Fp) or (N,win)
        flags = parse2inputs(flags,N,fp);

    case {3,4}
        % (N,DEV,'dev') or ('minorder',Fo,Dev) or ('minorder',Fo,Dev,'kaiser')
        flags = parse3inputs(flags,N,fp,varargin{:});
end

%------------------------------------------------------------------------
function bol = isvalidscalar(a)

bol = true;

if length(a) > 1 || ~isnumeric(a) || isnan(a) || isinf(a) || isempty(a)
    bol = false;
end
%------------------------------------------------------------------------
function [flags,flagschanged] = checklaststr(flags,str)

flagschanged = 0;

stringOpts = {'low','high','minphase','nonnegative'};
lpindx = find(strncmpi(str,stringOpts,length(str)));
if ~isempty(lpindx) && lpindx == 1
    flagschanged = 1;
elseif ~isempty(lpindx) && lpindx == 2
    flags.lowpassflag = 0;
    flagschanged = 1;
elseif ~isempty(lpindx) && lpindx == 3
    flags.minphase = 1;
    flagschanged = 1;
elseif ~isempty(lpindx) && lpindx == 4
    flags.zerophase = 1;
    flagschanged = 1;
end
function flags = parse2inputs(flags,N,fp)

if isnumeric(N) && length(fp) == 1 && ~flags.minphase && N < 6
    error(message('dsp:firhalfband:FilterErr2e'));
end
if dsp.util.isChar(N) && strncmpi(char(N),'minorder',length(char(N)))
    error(message('dsp:firhalfband:FilterErr1'));
end
checkvalidorder(N,flags);

if length(fp) > 1
    % Window specified, make sure it is of the right length
    if length(fp) ~= N + 1
        error(message('dsp:firhalfband:InvalidDimensions'));
    end
    flags.winflag = 1;
else
    if ~isvalidscalar(fp)
        error(message('dsp:firhalfband:FilterErr3'));
    end

    % Make sure fp is in the right range
    if fp <= 0 || fp >= 0.5
        error(message('dsp:firhalfband:InvalidRange'));
    end
end

%------------------------------------------------------------------------
function flags = parse3inputs(flags,N,fp,varargin)


if dsp.util.isChar(N)
    if strncmpi(char(N),'minorder',length(char(N)))
        flags.minordflag = 1;

        if ~isvalidscalar(varargin{1})
            error(message('dsp:firhalfband:FilterErr4'));
        end

        if nargin > 4 && dsp.util.isChar(varargin{2})
            % If a peak ripple is wanted, check if you want a Kaiser window design
            if strncmpi(char(varargin{2}),'kaiser', length(char(varargin{2})))
                flags.kaiserflag = 1;
                flags.winflag = 1;
            end
        end
    else
        error(message('dsp:firhalfband:FilterErr2a'));
    end

else
    checkvalidorder(N,flags);
    if ~dsp.util.isChar(varargin{1})
        error(message('dsp:firhalfband:FilterErr7'));
    end
    if ~strncmpi(char(varargin{1}),'dev', length(char(varargin{1})))
        error(message('dsp:firhalfband:FilterErr8'));
    end
    flags.firceqripflag = 1;
end

if ~isvalidscalar(fp)
    error(message('dsp:firhalfband:FilterErr3'));
end

% Make sure fp is in the right range
if fp <= 0 || fp >= 0.5
    error(message('dsp:firhalfband:InvalidRange'));
end

%------------------------------------------------------------------------
function checkvalidorder(N,flags)

if ~isvalidscalar(N)
    error(message('dsp:firhalfband:FilterErr2a'));
end
if ~flags.minphase && rem(N,2)
    error(message('dsp:firhalfband:FilterErr2b'));
elseif flags.minphase && ~rem(N,2)
    error(message('dsp:firhalfband:FilterErr2c'));
end
if flags.zerophase && rem(N,2)
    error(message('dsp:firhalfband:FilterErr2d'));
end

% [EOF]






