## Copyright (C) 2014   Robert Medina   <m05rob@aol.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/>.


## -*- texinfo -*-
## @deftypefn {taylorwin.m} {[@var{w}] =} taylorwin(@var{m})
## @deftypefnx {taylorwin.m} {[@var{w}] =} taylorwin(@var{m}, @var{n_bar})
## @deftypefnx {taylorwin.m} {[@var{w}] =} taylorwin(@var{m}, @var{n_bar}, @var{SLL})
##
##  m is the number of coefficients
##
##  n_bar is the number of nearly constant level sidelobes adjecent
##  to the mainlobe. Defaults to 8 if not provided.
##
##  SLL is the peak sidelobe level (in dB) relative to the mainlobe peak.
##  SLL must be a negative number. Defaults to -60dB if not provided.
##
##  length(w) == m
##
##  Creates a Taylor window coefficents
##
##  References:
##  This script was based on the equations that
##  are defined in the book titled Spotlight Synthetic
##  Aperture Radar Signal Processing Algorithms, By
##  Walter G. Carrara, Ron S. Goodman, and Ronald M. Majewski
##  ISBN 0-89006-728-7 From Pg 512 - 513.
## @end deftypefn

function out = taylorwin (m, n_bar, sll)

  ##Error Checking
  switch (nargin)
    case 1
      sll = -60;
      n_bar = 8;
    case 2
      if ( (n_bar <= 0) || (fix (n_bar) != n_bar) )
        error ("taylorwin: n_bar must be a postive integer greater than zero.")
      endif
      ##Set sll to default value
      sll = -60;
    case 3
      if(sll > 0)
        error ("taylorwin: SLL must be a negative value in dB");
      endif
      if ( (n_bar <= 0) || (fix (n_bar) != n_bar) )
        error ("taylorwin: n_bar must be a postive integer greater than zero.")
      endif
    otherwise
      print_usage();
  endswitch

  ##Always check m
  if ( (fix (m) != m) || (m <= 0) || ! (isscalar (m)) )
    error ("taylorwin: N must be a positive integer.")
  endif

  ##varaibles
  out = zeros(1, m);
  fq = zeros(1, n_bar - 1);

  ##temp variables
  temp = 0; 
  temp_sum = 0;

  ##Find the B Constant
  b = b_constant (sll);

  ##Find the A Constant
  a = a_constant (b);

  ##A squared
  a_sq = a * a;

  ##Get the value for Sigma^2
  sigma_sq = sigma (n_bar, a_sq);

  ##Find the values for fq
  for index = 1:n_bar - 1
    fq(index) = f_constant (sigma_sq, a_sq, n_bar, index, a);
  endfor

  ##Caculate Taylor coeffs
  ##Note out_index = m_index
  for out_index = 1:1:m

    ##adjust the out index var to start at 0, but matlab/octave arraies start
    ##at 1 so out_index is left at start value 1
    n = out_index - 1;
    temp_sum = 0;

    ##Find the sum of fq * Cos(arg)
    ##Note q = inner_index
    for inner_index = 1:n_bar - 1
      temp = fq(inner_index) * cos_arg(m, inner_index, n);
      temp_sum = temp + temp_sum;
    endfor

    temp_sum = temp_sum * 2;
    out(out_index) = 1 + temp_sum;

  endfor

endfunction

##  Set up the a constants for the window
##  The equation for the a constant is below
##  NOTE: log == ln
##       ln( b + sqrt ( b ^ 2 - 1 ) )
##  a = ----------------------------
##                  pi
##  @param b the b constant
##  @return returns a the constant
function a = a_constant (b)

  a = b * b;
  a = a - 1;
  a = sqrt(a);
  a = a + b;
  a = log(a);
  a = a / pi;

endfunction

##  Set up the b Constant
##  The equation for the b constant is below
##  b = 10 ^ ( - SLL / 20 )
##  @param SLL is the peak sidelobe level (in dB)
##     relative to the mainlobe peak
##  @return b the constant
function b = b_constant (sll)

  b = sll / 20;
  b = b * -1;
  b = 10^b;

endfunction

##  Set up sigma^2
##  The equation for sigma^e is below
##                         n_bar ^ 2
##   sigma^2 =  -------------------------------
##               A ^ 2 + ( n_bar - 0.5) ^ 2
##  @param n_bar the number of nearly constant level sidelobes adjecent
##     to the mainlobe
##  @param A_sq the value of the constant A^2
##  @return sigma_sq the constant
function sigma_sq = sigma (n_bar, a_sq)

  sigma_sq = n_bar * n_bar;
  temp = n_bar - 0.5;
  temp = temp * temp;
  temp = temp + a_sq;
  sigma_sq = sigma_sq / temp;

endfunction


## Find the value for the fq constant
##  @param sigma_sq
##  @param A_sq
##  @param n_bar
##  @param m the coeff number
##  @param A
##  @return fq
function fq = f_constant (sigma_sq, a_sq, n_bar, m, a)

  fq = top_multi (a, sigma_sq, m, n_bar, a_sq);
  fq = fq / bottom_multi (m, n_bar);

endfunction


##  Find the value for the numerator of fq,
##  Note: m > 1
##  The equation is below
##                        n_bar - 1
##                      ------------         m ^ 2 / sigma_sq
##  out =  (-1)^(m+1)   |          |   1 - -----------------------------
##                      |          |         A_sq + ( index - 0.5 ) ^ 2
##                        index = 1
##
##
function out = top_multi (a, sigma_sq, m, n_bar, a_sq)

  #varaible
  m_sq = m * m;

  top_frac = m_sq / sigma_sq;
  temp = 0;
  out = 1;

  for index = 1:1:n_bar - 1

    bottom = index - 0.5;
    bottom = bottom * bottom;
    bottom = bottom + a_sq;
    temp = top_frac / bottom;
    temp = 1 - temp;
    out = out * temp;

  endfor

  if ( mod( (m + 1), 2 ) > 0 )
    out = out * -1;
  endif

endfunction


##  Find the value for the numerator of fq,
##  Note: m > 1
##  The equation is below
##             n_bar - 1
##           ------------           m ^ 2
##  out = 2  |          |   1 -  -----------
##           |          |         index ^ 2
##            index = 1
##            index =/= m
##
function out = bottom_multi (m, n_bar)

  m_sq = m * m;
  out = 1;
  temp = 0;

  for index = 1:1:n_bar - 1

    if (index == m)
      continue;
    endif

    temp = index * index;
    temp = m_sq / temp;
    temp = 1 - temp;
    out = out * temp;

  endfor

  out = out * 2;

endfunction


##  Find the cosine function value of the Taylor window.
##  The equation for the cosine function for the taylor window is define below
##
##       2 * pi * (n - 0.5 * N + 0.5)
##  cos( ---------------------------- )  = out
##                   N
##
##  @param n_size the total number of coefficents
##  @param m the value for fq (Or index value)
##  @param n the number of the coefficent that is being currently caculated
##  @return out the value of the cosine function defined above
function out = cos_arg (n_size, m, n)

  out = n_size * 0.5;
  out = n - out;
  out = out + 0.5;
  out = out * 2 * pi * m;
  out = out / n_size;
  out = cos(out);

endfunction

#Catch no Args
#taylorwin();

#Catch Non pos int for size
#taylorwin(-1);
#taylorwin(-1,8);
#taylorwin(-1,8,-60);

#Catch Non int for size
#taylorwin(0.5);
#taylorwin(0.5,8);
#taylorwin(0.5,8,-60);

#Catch Non Pos int for N_Bar
#taylorwin(171,-1);
#taylorwin(171,0.5);

#Catch Non Neg SLL
#taylorwin(171,8,60);
