function [n, fo, ao, w] = firpmord(f, a, dev, varargin)
// Parks-McClennan optimal FIR filter order estimation
//
//
// Calling sequence
// [n,fo,ao,w] = firpmord(f,a,dev)
// [n,fo,ao,w] = firpmord(f,a,dev,fs)
//
// 
// Parameters
// f: double - positive - vector
//      Frequency band edges (between 0 and Fs/2).
//      Length of f is two less than the length of a.
// a: double - positive - vector
//      Desired amplitudes at the frequency bands specified by f.
// dev: double - positive - vector
//      Maximum allowable deviations.
//      Maximum acceptable deviations or ripples between the frequency response
//      and the desired amplitudes in the frequency bands specified by f. Must have 
//      the same length as a.
// n: int - scalar
//      Filter order
// fo: double - positive - vector
//      Frequency vector
// ao: double - positive - vector
//      Amplitude vector
// w: double - vector
//      Weights
//
// 
// Examples
// [1] A low-pass filter
//          f = [1500 2000];        // frequency edges for bands
//          a = [1 0];              // desired amplitude for each band
//          dev = [0.01 0.1];       // Acceptable deviation for each band
//          fs = 8000;              // Sampling frequency
//          [n,fo,ao,w] = firpmord(f,a,dev,fs);
//
// [2] A bandstop filter
//          f = [1000 1800 2400 3000];
//          a = [1 0 0.5];
//          dev = [0.01 0.1 0.03];
//          fs = 8000;
//          [n,fo,ao,w] = firpmord(f,a,dev,fs);
//
//
// References
// [1] Rabiner, Lawrence R., and Bernard Gold. "Theory and application of 
//     digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc., 
//     1975. 777 p. 156-7 (1975).
// [2] Rabiner, Lawrence R., and Otto Herrmann. "The predictability of certain 
//     optimum finite-impulse-response digital filters." Circuit Theory, 
//     IEEE Transactions on 20.4 (1973): 401-408.
//
// Authors
// Ayush Baid
//
//
// See Also
// buttord | cheb1ord | cheb2ord | ellipord | firpm | kaiserord

    [numOutArgs,numInArgs] = argn(0);
    
    
    // ********************
    // Checking number of arguments
    // ********************
    
    if numInArgs~=3 & numInArgs~=4 then
        msg = "firpmord: Wrong number of input argument; 3-4 expected";
        error(77,msg);
    end
    
    if numOutArgs~=4 then
        msg = "firpmord: Wrong number of output argument; 4 expected";
        error(78,msg);
    end
    
    // ********************
    // Parsing input args
    // ********************

    // Parsing fs
    fs = 2; // default
    if length(varargin)==1 then
        fs = varargin(1);
        if length(fs)~=1 then
            msg = "firpmord: Wrong type for argument #4 (fs): Positive real scalar expected";
            error(53,msg);
        end
        if fs<=0 then
            msg = "firpmord: Wrong type for argument #4 (fs): Positive real scalar expected";
            error(53,msg);
        end 
        if type(fs)~=1 & type(fs)~=8 then
            msg = "firpmord: Wrong type for argument #4 (fs): Positive real scalar expected";
            error(53,msg);
        end
    end
    
    // Checks on f
    if ~isvector(f) | (type(f)~=1 & type(f)~=8) then
        msg = "firpmord: Wrong type for argument #1 (f): Vector of positive reals expected";
        error(53,msg);
    end
    
    if ~(and(f>=0) & and(f<=fs/2)) then
        msg = "firpmord: Wrong value for argument #1 (f): Values must be between 0 and fs/2";
        error(116,msg);
    end
    
    
    // Check on a
    if ~isvector(a) | (type(a)~=1 & type(a)~=8) then
        msg = "firpmord: Wrong type for argument #2 (a): Vector of positive reals expected";
        error(53,msg);
    end
    if ~and(a>=0) then
        msg = "firpmord: Wrong value for argument #2 (a): Values must be positive";
        error(116,msg);
    end
    
    if length(f)~=2*length(a)-2 then
        msg = "firpmord: Wrong type for arguments #1(f) and #2 (a): Length of f must be two less than twice the length of a ";
        error(53,msg);
    end
    
    
    // Check on dev
    if ~isvector(dev) | (type(dev)~=1 & type(dev)~=8) then
        msg = "firpmord: Wrong type for argument #3 (dev): Vector of positive reals expected";
        error(53,msg);
    end
    if ~and(dev>0) then
        msg = "firpmord: Wrong value for argument #3 (dev): Values must be positive";
        error(116,msg);
    end
    if length(dev)~=length(a) then
        msg = "firpmord: Wrong type for arguments #2(a) and #3 (dev): Length of a and dev must be equal";
        error(53,msg);
    end
    
    // ********************
    // Some preprocessing
    // ********************
    
    // Turn every vector into a column vector
    f = f(:);
    a = a(:);
    dev = dev(:);
    
    // Normalizing frequencies
    f = f./fs;
    
    // Get deviation relative to the amplutudes
    is_zero = a==0; 
    dev = dev./(is_zero+a); // no scaling req. when desired amplitude is 0
    
    
    num_bands = size(a,1);
    
    // Dividing frequency band edges into 2 vectors, f1 and f2, denoting 
    // passband and stopband edges respectively.
    f1 = f(1:2:$-1);
    f2 = f(2:2:$);
    
    
    // ********************'
    // Calculations for filter order
    // ********************
    
    // Note: Amplitudes don't matter for order as they can be adjusted by 
    //       scaling and linear shifting    
    
    if num_bands==2 then
        // Simple low-pass or high-pass filter, use single_transition_order_estimation
        
        L = single_trans_order_est(f1(1), f2(1), dev(1), dev(2));
    else
        // We have a bandpass filter, which will be considered to be composed
        // of a cascade of simple high-pass/low-pass filters
        
        // The first filter is considered high-pass (if it is low pass in reality,
        // the filter just has to be negated; filter order does not change)
        
        // Loop over different simple filters and select the highest required length
        L = 0
        for i=2:num_bands-1
            L1 = single_trans_order_est(f1(i-1), f2(i-1), dev(i), dev(i-1));
            L2 = single_trans_order_est(f1(i), f2(i), dev(i), dev(i+1));
            L = max([L; L1; L2]);
        end
    end
    
    // filt. order = L-1
    n = ceil(L) - 1;
    
    // frequency and respective amplitudes
    fo = [0;2*f;1];
    ao = zeros(size(fo,1),1);
    ao(1:2:$) = a;
    ao(2:2:$) = a;
    
    // weights
    w = ones(size(dev,1),1)*max(dev)./dev;            
    
endfunction

function L = single_trans_order_est(freq1, freq2, delta1, delta2)
// Calculates the filter order for a single transition band filter ( simple 
// low-pass or simple high-pass filter)
//
//
// Parameters (assuming a low-pass filter; notations change for high pass filter)
//      freq1: passband cutoff frequency (normalized)
//      freq2: stopband cutoff frequency (normalized)
//      delta1: passband ripple (max. allowed)
//      delta2: stopband attenuation (not in dB)
//      L: filter length; filter order = L-1
//
// Note: will not work well when transition near f = 0 or 0.5
//
// References
// [1] Rabiner, Lawrence R., and Bernard Gold. "Theory and application of 
//     digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc., 
//     1975. 777 p. 156-7 (1975).

    // Creating a matrix for consice representation of coeffs a_i used in Eq. 3.142
    // and b_i in Eq. 3.143 in [1]
    A = [-4.278e-1, -4.761e-1; -5.941e-1, 7.114e-2; -2.66e-3, 5.309e-3];
    B = [11.01217, 0.51244];

    log_delta1 = log10(delta1);
    log_delta2 = log10(delta2);

    // Evaluating Eq. 3.142 (Ref. [1])
    D = [1, log_delta1, log_delta1.^2] * A * [1; log_delta2];

    // Evaluating Eq. 3.143 (Ref. [1])
    f = B * [1; log_delta1 - log_delta2];

    // Evaluating Eq. 3.145 (Ref. [1])
    del_freq = abs(freq1 - freq2);
    L = D./del_freq - f.*del_freq + 1;

endfunction