## 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{N})
## @deftypefnx {taylorwin.m} {[@var{w}] =} taylorwin(@var{N}, @var{n_bar})
## @deftypefnx {taylorwin.m} {[@var{w}] =} taylorwin(@var{N}, @var{n_bar}, @var{SLL})
##
##  N 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) == N
##  
##  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(N, n_bar, SLL)

   #Error Checking
   if ( nargout > 1 )
      print_usage;
   endif
   
   switch(nargin)
      case 1
         SLL = -60;
         n_bar = 8;
      case 2
         if ( (n_bar <= 0) || (fix(n_bar) != n_bar) )
            error ("n_bar must be a postive integer greater than zero.")
         endif
         #Set SLL to default value
         SLL = -60;
      case 3
         if(SLL > 0)
            error ("SLL must be a negative value in dB");
         endif
         if ( (n_bar <= 0) || (fix(n_bar) != n_bar) )
            error ("n_bar must be a postive integer greater than zero.")
         endif
      otherwise
        print_usage;
   endswitch
   
   if (fix(N) != N)
      N = round(N);
      warning("N must be an integer, so rounding N to the nearest integer.")
   endif

   #varaibles
   out = zeros(1,N);
   Fm = 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 Fm
   for index = 1:1:n_bar - 1

      #           F_Constant(sigma_sq, A_sq, n_bar, m, A)
      Fm(index) = F_Constant(sigma_sq, A_sq, n_bar, index, A);

   endfor

 
   #Caculate Taylor coeffs
   #Note out_index = n
   for out_index = 1:1:N

      #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 Fm * Cos(arg)
      # Cos_arg(N, m, n)
      #Note m = inner_index
      for inner_index = 1:1:n_bar - 1
         temp = Fm(inner_index) * Cos_arg(N, inner_index, n);
         temp_sum = temp + temp_sum;
      endfor

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

   endfor

end

#  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;

end

#  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;

end

#  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;

end


#Find the value for the Fm constant
#  @param sigma_sq
#  @param A_sq
#  @param n_bar
#  @param m the coeff number
#  @param A
#  @return Fm
function [Fm] = F_Constant(sigma_sq, A_sq, n_bar, m, A)

   #    Top_Multi(A, sigma_sq, m, n_bar, A_sq)
   Fm = Top_Multi(A, sigma_sq, m, n_bar, A_sq);
   Fm = Fm / Bottom_Multi(m, n_bar);

end


#  Find the value for the numerator of Fm, 
#  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

end



#  Find the value for the numerator of Fm, 
#  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;

end


#  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 the total number of coefficents
#  @param m the value for Fm (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, m, n)

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

end
