## Copyright (C) 2013 James-Adam Renquinha Henri <jarhmander@gmail.com>
##
## 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/>.

## Generate poles, zero and gain of an analog prototype Legendre-Papoulis filter
## (also known as Optimal "L" filter).
##
## Compared to Butterworth filters, they have a much sharper rolloff at the 
## expense of a somewhat less flat passband; but it has monotonic frequency
## response unlike the Chebyshev type I.
##
## [z,p,k] = legendap(N)
##    analog prototype filter of order N

function [z,p,k] = legendap(N)
    if N < 1 || N > 50, error('Invalid order'), endif
    z = [];
    k = ceil(N/2) - 1;  % k here is upper bound of sum
    p = 0;
    isNodd = mod(N,2) == 1;
    ai=[];
    
    ## Define ai according to parity of N and possibly k
    if isNodd
        pden = 1/(sqrt(2)*(k + 1));
        ai = @(i) (pden*(2*i+1));
    else
        pden = 1/sqrt((k+1)*(k+2));
        Kmod2 = mod(k,2);
        ai = @(i) (pden*((i*2+1)*(mod(i,2) == Kmod2)));
    endif
        
    ## Accumulate polys
    for i=0:k
        p = polyadd(p, ai(i) * legendrepoly(i));
    endfor
    
    ## Square poly
    p = conv(p,p);
    if ~isNodd
        p = conv([1 1], p);
    endif
    p = polyint(p);
    
    ## This is to evaluate the analytical integral, (x -> w^2 -1)
    p = polysub(polysubst(p, [2 0 -1]), polyval(p, -1));
    
    ## w^2 -> -s^2
    p = upsample(polysubst(downsample(p, 2), [-1 0]),2);
    p = p(1:end-1);
    p(end) = 1;
    
    ## The coefficients of the "L" polynomials are integers, so the rounding here
    ## will fix small numerical inaccuracies.
    p = round(p);
    
    ## Get the roots and only keep those with negative real parts
    p = roots(p);
    p = p(real(p) < 0);    
    
    ## put the correct dc gain
    k = prod(abs(p)); 
endfunction

function y = polyadd(w,x)
    m = min(length(w),length(x));
    if length(w) > length(x), y=w; else y=x; end
    y(end-m+1:end) = w(end-m+1:end) + x(end-m+1:end);
endfunction

function y = polysub(w,x)
    y = polyadd(w,-x);
endfunction

function y = polysubst(p,v)
    ip = 1;
    y = 0;
    for c = p(end:-1:1)
        y = polyadd(y, c*ip);
        ip = conv(ip,v);
    endfor
endfunction
