(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)


(* my Mathematica package, contains collection of functions for different purposes.
See test_mma.nb for documenations and how to use each function.
By Nasser M. Abbasi

January 9, 2017   added  checkForSingularity
Nov 29. 2014      started
*)

BeginPackage["nma`"]

Unprotect@@Names["nma`*"];
ClearAll@@Names["nma`*"];

getStateGainVector::usage="getStateGainVector[A,B,designEigenvalues,debugFlag]

Takes state space A,B, as in x'=Ax+Bu and list of design eigenvalues
and returns state feedback gain vector k, such that A+B*k have the desired eigenvalues.

If debugFlag is True, then design steps are printed.

Example:

a={{0,0,1,0},{0,0,0,1},{-2,1,0,0},{1,-1,0,0}};
b={{0},{0},{1},{0}};
desiredEigenValues={-1+I,-1+I,-1-I,-1-I};
k = getStateGainVector[a,b,desiredEigenValues,True]
   {-5,2,-4,-4}

"


getControllableForm::usage="getControllableForm[A,B]

Takes A,B state space matrices and returns the controllable canonical form.
Example

a={{0,0,1,0},{0,0,0,1},{-2,-1,0,0},{1,-1,0,0}};
b={{0},{0},{1},{0}};
{A0,B0}=getControllableForm[a,b]

"


generateLTVObservabilityMatrix::usage="generateLTVObservabilityMatrix[A,C,t]

Generate the Linear time varying observability matrix M defined as 

M0=C
M(k+1)= M(k) A(t) + d/dt M(k)  for k=0..n-2

where n is the size of the square A matrix.

If the rank of the matrix is n at t>t0, where n is the size of the square matrix A,
then the system is observable at t0.

Example:

a={{0,1,t},{t^2,-t,Exp[t]},{1,-2,1}};
c={{1,t,Exp[-t]},{1,1,2}};
(r=generateLTVObservabilityMatrix[a,c,t])//MatrixForm
MatrixRank[r/.t\[Rule]1]

"


generateLTVControllabilityMatrix::usage="generateLTVControllabilityMatrix[A,B,t]

Generate the Linear time varying controllability matrix M defined as 

M0 = B(t)
M(k+1) = -A(t) M(k) + d/dt M(k)   for k=0..n-2

If the rank of the matrix is n at t>t0, where n is the size of the square matrix A,
then the system is controllable at t0.

Example

a={{Cos[t],t},{Exp[t],2}};
b={{1},{t}};
r=generateLTVControllabilitymatrix[a,b,t]

 {{1, -t^2 - Cos[t]}, {t, 1 - E^t - 2*t}}

Det[r]

   1 - E^t - 2*t + t^3 + t*Cos[t]

Now evaluate the about for any t>t0 where t0 is where you want to check contrability at.

If the matrix r is not square, then use MatrixRank[r/.t\[Rule]1] and check for rank=n
"

displayHurwitz::usage="displayHurwitz[p_List,s_] 

p is list of polynomials in s.
Generates Hurwitz stability Matrix and displays it.

Example:

displayHurwitz[{11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5}, s]
"


hurwitz::usage="hurwitz[p_,s_]

Take one polynomial in s and reurns the Hurwitz stability matrix.

Example:
hurwitz[11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5, s]; 
Grid[%, Frame -> All]
"


leadingDet::usage="leadingDet[h_]

Take a Hurwitz matrix returned by hurwitz[polynomial,s] and returns the leading
minors determinants. They must all be postive for stable polynomial.

Example

hurwitz[11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5, s]
leadingDet[%]
    {9,6,-108,-196,-196}
"


kharitonovPoly::usage="kharitonovPoly[w_List,s_]

Take list of intervals and returns the four Kharitonov Polynomials in s.
Example, 

(nma`kharitonovPoly[{{1,2},{9,10}},s])//TableForm
\!\(\*TagBox[TagBox[GridBox[{
{
RowBox[{\"1\", \"+\", 
RowBox[{\"9\", \" \", \"s\"}]}]},
{
RowBox[{\"2\", \"+\", 
RowBox[{\"10\", \" \", \"s\"}]}]},
{
RowBox[{\"2\", \"+\", 
RowBox[{\"9\", \" \", \"s\"}]}]},
{
RowBox[{\"1\", \"+\", 
RowBox[{\"10\", \" \", \"s\"}]}]}
},\nGridBoxAlignment->{\"Columns\" -> {{Left}}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {{Baseline}}, \"RowsIndexed\" -> {}},\nGridBoxSpacings->{\"Columns\" -> {Offset[0.27999999999999997`], {Offset[0.5599999999999999]}, Offset[0.27999999999999997`]}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, \"RowsIndexed\" -> {}}],
Column],
Function[BoxForm`e$, TableForm[BoxForm`e$]]]\)
Notice that the order for input is like this: limits of a0 followed by limits of a1 followed by limits of a2 and so on.
Hence to find Kharitonov polynomials for  [1,2]+[9,10]s  we enter it as shown above.

Example 2: For polynomial [1,2]+[9,10]s+[-1,3]s^2 the call is

(nma`kharitonovPoly[{{1,2},{9,10},{-1,3}},s])//TableForm

\!\(\*TagBox[TagBox[GridBox[{
{
RowBox[{\"1\", \"+\", 
RowBox[{\"9\", \" \", \"s\"}], \"+\", 
RowBox[{\"3\", \" \", SuperscriptBox[\"s\", \"2\"]}]}]},
{
RowBox[{\"2\", \"+\", 
RowBox[{\"10\", \" \", \"s\"}], \"-\", SuperscriptBox[\"s\", \"2\"]}]},
{
RowBox[{\"2\", \"+\", 
RowBox[{\"9\", \" \", \"s\"}], \"-\", SuperscriptBox[\"s\", \"2\"]}]},
{
RowBox[{\"1\", \"+\", 
RowBox[{\"10\", \" \", \"s\"}], \"+\", 
RowBox[{\"3\", \" \", SuperscriptBox[\"s\", \"2\"]}]}]}
},\nGridBoxAlignment->{\"Columns\" -> {{Left}}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {{Baseline}}, \"RowsIndexed\" -> {}},\nGridBoxSpacings->{\"Columns\" -> {Offset[0.27999999999999997`], {Offset[0.5599999999999999]}, Offset[0.27999999999999997`]}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, \"RowsIndexed\" -> {}}],
Column],
Function[BoxForm`e$, TableForm[BoxForm`e$]]]\)
"


polarForm::usage="polarForm  wrapper to make polar form for complex number
written by Xerxes. see

http://mathematica.stackexchange.com/questions/20570/wrapper-for-inexact-numeric-complex-numbers-that-maintains-polar-form

Sqrt[5] E^(I ArcTan[2]) // N // polarForm
(* 2.23607E^(I1.10715) *)
"

(*--------------------------------------------*)
(*  writtenby Hypnotoad                *)
(* http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245*)
(*--------------------------------------------*)      
getMaxPadding::usage="getMaxPadding[p_List] 

finds maximum image padding for 
a list of plots. 

writtenby Hypnotoad from stackoverflow
http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245

p1=Plot[Sin[x],{x,-Pi,Pi},ImagePadding->30]
p2=Plot[Cos[x],{x,-Pi,Pi},ImagePadding->20]
Grid[{{p2, p1}} /. Graphics[x__] :> 
    Graphics[x, ImagePadding -> getMaxPadding[{p1, p2}]]]"

(*--------------------------------------------*)
(*  writtenby Hypnotoad                *)
(* http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245*)
(*--------------------------------------------*)          
snip::usage="snip[pos_] 

writtenby Hypnotoad
http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245

The function snip[x] places the mark on the axes at relative 
position x (0 and 1 being the ends)

example:

data1 = {{1, 1.1}, {2, 1.5}, {3, 0.9}, {4, 2.3}, {5, 1.1}};
data2 = {{1, 1001.1}, {2, 1001.5}, {3, 1000.9}, {4, 1002.3}, {5, 
    1001.1}};
p1 = ListPlot[data1, PlotRange -> All, Joined -> True, Mesh -> Full, 
   PlotStyle -> Red, AxesStyle -> {None, snip[1]}, 
   PlotRangePadding -> None, ImagePadding -> 40];
p2 = ListPlot[data2, PlotRange -> All, Joined -> True, Mesh -> Full, 
   PlotStyle -> Blue, Axes -> {False, True}, 
   AxesStyle -> {None, snip[0]}, PlotRangePadding -> None, 
   ImagePadding -> 40];

Column[{p2, p1} /. 
  Graphics[x__] :> 
   Graphics[x, ImagePadding -> getMaxPadding[{p1, p2}], 
    ImageSize -> 400]]"
    
(*--------------------------------------------*)
(*  written by John Fultz, see                *)
(* http://mathematica.stackexchange.com/questions/1039/mouseappearance-and-cursor-problems*)
(*--------------------------------------------*)  
stripGraphMouseAppearance::usage="stripGraphMouseAppearance[x_Graph]
puts normal arrow looking mouse appearance on Grph plots. 

 written by John Fultz, see
see http://mathematica.stackexchange.com/questions/1039/mouseappearance-and-cursor-problems
 example:
 stripGraphMouseAppearance[
  TreeGraph[{1 -> 2, 1 -> 3, 2 -> 4, 2 -> 5}, 
  VertexLabels -> {5 -> Placed[\"five\", Tooltip]}, ImagePadding -> 20]
 ]"
 
(*----------------*)
(*                *)
(*----------------*)
matrix2DredBlackPositions::usage="matrix2DredBlackPositions[nRow,nCol,opt] returns {r,b}, where 
r is a list of indices in a matrix of dimensions[nRow,nCol] where the red balls are, and 
b is list of indices where black balls are. These are the red/black positions to use 
for Gauss-Seidel red, black algorithm.

opt: internal->True, is the default, which says to apply red/black only to the 
internal nodes, and not to the boundaries.  

Example, put a '1' where a red should be, and put a '-1' where a black should be

AnRow = 3;
nCol  = 3;
a     = Table[0,{i,nRow},{j,nCol}];
{r,b}=nma`matrix2DredBlackPositions[nRow,nCol];
(a[[Sequence@@#]]=1)&/@r;
(a[[Sequence@@#]]=-1)&/@b;
MatrixForm[a]

"

(*----------------*)
(*                *)
(*----------------*)
matrix2DLexOrdering::usage="matrix2DLexOrdering[mat] returns {A,B}, where A is matrix 
of the same as as mat, which contains in each entry, the coordinates of each cell in the 
matrix in Lexicographical order and B contains the sequentional Lexicographical order

Example

a=Table[0,{i,3},{j,3}];
nma`matrix2DLexOrdering[a]
Out[310]= {{{1,3},{2,3},{3,3}},{{1,2},{2,2},{3,2}},{{1,1},{2,1},{3,1}}}

"

(*----------------*)
(*                *)
(*----------------*)
centerImageData::usage="centerImageData[im] centers the image by multiplying each pixel by (-1)^(i+j) 
where i,j are the row number and column number of the pixel in normal matrix meaning of row and column. 
This is useful to do before taking the FFT of the image.
"

(*----------------*)
(*                *)
(*----------------*)
diagonalBlockMatrix::usage="diagonalBlockMatrix[block,n] gives a diagonal square matrix 
with the block along the diagonal repeated n times

example:

In[14]:= block={{1,2},{3,4}}
Out[14]= {{1,2},{3,4}}

Normal@nma`diagonalBlockMatrix[block,4] 
Out[15]= {{1,2,0,0,0,0,0,0},
          {3,4,0,0,0,0,0,0},
          {0,0,1,2,0,0,0,0},
          {0,0,3,4,0,0,0,0},
          {0,0,0,0,1,2,0,0},
          {0,0,0,0,3,4,0,0},
          {0,0,0,0,0,0,1,2},
          {0,0,0,0,0,0,3,4}}

"

(*----------------*)
(*                *)
(*----------------*)
poissonMatrixOnSquareDirichlet::usage="poissonMatrixOnSquareDirichlet[n,m] gives the 
discretization A matrix for \!\(\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\)u+\!\(\*SubscriptBox[\(\[PartialD]\), \(y, y\)]\)u 
using Laplacian stensil of size m (currently only m=5 is supported. n is the number of grid points along one 
dimension of the square. Hence, for n=3, this means the square has 9 internal grid points. These are 
the internal grid points only, as this is Dirichlet boundary conditions problem

Example

nma`poissonMatrixOnSquareDirichlet[3,5]//MatrixForm

"

(*----------------*)
(*                *)
(*----------------*)
diagonalDominantMatrixQ::usage="diagonalDominantMatrixQ[A] gives True if matrix is diagonal dominant matrix, and False otherwise. The following are available options

strict->{True,False}

strict->False, is the default. The check is made using >= for dominance of diagonal. 

stict->True, the check is made using > for dominance of diagonal

Examples:

<<nma`
mat={{4,4},{4,4}} 

nma`diagonalDominantMatrixQ[mat,strict->False]
Out[27]= True

nma`diagonalDominantMatrixQ[mat,strict->True]
Out[31]= False

nma`diagonalDominantMatrixQ[mat]
Out[38]= True

"

(*----------------*)
(*                *)
(*----------------*)
numItse::usage="numItse[num,w,d,fontsize] formats a number in scientific form. 
num is the number to format
w is the total width of the field including the decimal point
d is width of field to right of decimal point
fontsize is the font size to use
 
EXAMPLE
nma`numItse[234.4,6,4,10,12]

\!\(\*
StyleBox[\"2.3440\",\nStripOnInput->False,\nFontSize->14]\)\!\(\*
StyleBox[\"e2\",\nStripOnInput->False,\nFontSize->14]\)    
 "
(*----------------*)
(*                *)
(*----------------*)
linspace::usage="generate linearly spaced list.
    linespace(x1, x2) generate list of 100 equally spaced points between x1 and x2.
 
    linspace(x1, x2, N) generate N equally spaced points between x1 and x2. if N < 2, it returns x2.
 "
(*----------------*)
(*                *)
(*----------------*)
str::usage="str[expr] removes $nnn from the expression. Used for printing from inside a function"
(*----------------*)
(*                *)
(*----------------*)
formatBinary::usage="formatBinary[num,nBits]
returns a string representation of the number num in base 2. nBits is the number of bits available to use for displaying the final result.

EXAMPLE
nma`formatBinary[4,4]
0100

nma`formatBinary[4,8]
0000 0100
"
(*----------------*)
(*                *)
(*----------------*)
toTwoComplement::usage="toTwoComplement[x,mp,nBits] returns {delta,xQuantized,code}
This function calculates the 2 complement value (in decimal) of the analog value x after quantization. mp is the maximum level allowed for the quantization. nBits is the number of bits allowed to store the final result in. The code is the two's complement value generated. To convert this value to binary reprentation, use the funtion formatBinary[] in this package.
"
(*----------------*)
(*                *)
(*----------------*)
toOnesComplement::usage="toOnesComplement[x,mp,nBits] returns {delta,xQuantized,code}
This function calculates the ones complement value (in decimal) of the analog value x after quantization. mp is the maximum level allowed for the quantization. nBits is the number of bits allowed to store the final result in. The code is the ones complement value generated. To convert this value to binary reprentation, use the funtion formatBinary[] in this package."
(*----------------*)
(*                *)
(*----------------*)
toSignMagnitude::usage="toSignMagnitude[x,mp,nBits] returns {delta,xQuantized,code}
This function calculates the sign magnitude value (in decimal) of the analog value x after quantization. mp is the maximum level allowed for the quantization. nBits is the number of bits allowed to store the final result in. The code is the sign magnitude value generated. To convert this value to binary reprentation, use the funtion formatBinary[] in this package.

Example

nBits=8;
mp=16;
x= 8.7;
       {delta,xQuantized,code}=toSignMagnitude[x,mp,nBits]

"

(*----------------*)
(*                *)
(*----------------*)
toOffsetBinary::usage="toOffsetBinary[x,mp,nBits] returns {delta,xQuantized,code}

This function calculates the offset binary value (in decimal) of the analog value x after quantization. mp is the maximum level allowed for the quantization. nBits is the number of bits allowed to store the final result in. The code is the offset binary value generated. To convert this value to binary reprentation, use the funtion formatBinary[] in this package.

Example

nBits=8;
mp=16;
x=-8.7;
         {delta,xQuantized,code}=toOffsetBinary[x,mp,nBits]

"

(*----------------*)
(*                *)
(*----------------*)
quantizationLevel::usage="quantizationLevel[x,mp,nBits]

Find the quantization level for an analog value x. mp is the maximum allowed level and nBits is the number of bits available for encoding x.

Example

nBits=8;
mp=16;
x=-6.5;

quantizationLevel[-8.7,16,nBits]

58
"

(*----------------*)
(*                *)
(*----------------*)
getPhaseOfPulseTrainFourierCoeff::usage="getPhaseOfPulseTrainFourierCoeff[delay,period,range,dutyCycle,nCoeff]

Example

    nma`getPhaseOfPulseTrainFourierCoeff[0,1,10,1,5]

           {{-4,0},{-3,0},{-2,0},{-1,0},{0,0},{1,0},{2,0},{3,0},{4,0}}

"

(*----------------*)
(*                *)
(*----------------*)
getMagnitudeOfPulseTrainFourierCoeff::usage="getMagnitudeOfPulseTrainFourierCoeff[t0,period,range,dutyCycle,nCoeff,h]

Example

nma`getMagnitudeOfPulseTrainFourierCoeff[0,1,10,1,5,3]

    {{-4,0},{-3,0},{-2,0},{-1,0},{0,3},{1,0},{2,0},{3,0},{4,0}}

"

(*----------------*)
(*                *)
(*----------------*)
makePulseTrain::usage="makes a pulse train. Returns List of Lines that can be used in a plot

makePulseTrain[delay_?(NumberQ[#]&&Positive[#]&),(* delay of a pulse from time=0*)
period_?(NumberQ[#]&&Positive[#]&),  (* the time between each pulse*)
range_?(NumberQ[#]&&Positive[#]&),   (* the time duration to use*)
tao_?(NumberQ[#]&&Positive[#]&),     (* width of each pulse*)
h_?(NumberQ[#]&&Positive[#]&)        (* height of the pulse*)
]

Returns List of Lines.   

EXAMPLE USAGE:

delay=.2;period=1;range=5.;tao=.5; h=1;

currentPulses=nma`makePulseTrain[delay,period,range,tao,h];

Plot[0,{x,-range,range},
       PlotRange\[Rule]{{-range,range},{0,h}},
       Epilog\[Rule]{Thin,Red,currentPulses}
    ]
"
(*----------------*)
(*                *)
(*----------------*)
fcPulseTrain::usage="fcPulseTrain[n,h,width,period]
returns the n'th fourier series complex coeff for a pulse train of height h, of width and period

EXAMPLE:

h=1;width=1;period=2;
 
Table[{n,nma`fcPulseTrain[n,h,width,period]},{n,-4,4}]//TableForm

\!\(\*TagBox[GridBox[{
{\(-4\), \(0\)},
{\(-3\), \(-\*FractionBox[\(1\), \(3\\\ \[Pi]\)]\)},
{\(-2\), \(0\)},
{\(-1\), FractionBox[\(1\), \(\[Pi]\)]},
{\(0\), FractionBox[\(1\), \(2\)]},
{\(1\), FractionBox[\(1\), \(\[Pi]\)]},
{\(2\), \(0\)},
{\(3\), \(-\*FractionBox[\(1\), \(3\\\ \[Pi]\)]\)},
{\(4\), \(0\)}
},\nGridBoxAlignment->{\"Columns\" -> {{Left}}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {{Baseline}}, \"RowsIndexed\" -> {}},\nGridBoxSpacings->{\"Columns\" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, \"RowsIndexed\" -> {}}],
Function[BoxForm`e$, TableForm[BoxForm`e$]]]\)
"

(*----------------*)
(*                *)
(*----------------*)
plotFourierTransform::usage="plotFourierTransform[ft,w,from,to,ymin,ymax,arrowheadsize]
Plot fourier transform, including pulses

EXAMPLE

x[t_]:=2 Cos[2 Pi 10 t]
Plot[x[t],{t,0, 2/10}]
FourierTransform[x[t],t,w,FourierParameters\[Rule]{-1,-1}];

nma`plotFourierTransform[%,t,0,2 period,-2,2,Small]

"


(*----------------*)
(*                *)
(*----------------*)
numIt::usage="Function takes a number and return it as string formatted

numIt[v_?(NumberQ[#]&),    (*the number to format*)
s1_?(IntegerQ[#]&& Positive[#] &), (* total number of digits to use, including the decimal digits*)
s2_?(IntegerQ[#]&& Positive[#] &),  (*number of digits to display to the right of the decimal point*)
fontSize_?(IntegerQ[#]&& Positive[#] &) (*font size to use for string formating*)

EXAMPLE USAGE:

nma`numIt[23.321,5,1,14]

23.3

"

(*----------------*)
(*                *)
(*----------------*)
butterToLowPass::"usage"="butterToLowPass[Hs,\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(c\)]\),variable];
convert Butterworth normalized low pass analog filter Hs to low pass analog filter with cutoff frequency \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(c\)]\). Variable is the symbolic variable of the polynomial.
"

(*----------------*)
(*                *)
(*----------------*)
butterToHighPass::"usage"="butterToHighPass[Hs,\[CapitalOmega],variable];
convert Butterworth normalized low pass analog filter Hs to high pass analog filter with corner frequency \[CapitalOmega]. Variable is the symbolic variable of the polynomial.
"

(*----------------*)
(*                *)
(*----------------*)
butterToBandPass::"usage"="butterToBandPass[Hs,\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\),variable];
convert Butterworth normalized low pass analog filter Hs to band pass analog filter with low corner frequency \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\) and high corner frequency \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\). Variable is the symbolic variable of the polynomial.
"

(*----------------*)
(*                *)
(*----------------*)
butterToBandStop::"usage"="butterToBandStop[Hs,\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\),variable];
convert Butterworth normalized low pass analog filter Hs to band stop analog filter with low corner frequency \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\) and high corner frequency \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\). Variable is the symbolic variable of the polynomial.
"

(*----------------*)
(*                *)
(*----------------*)

getButterworthPolynomial::"usage"="getButterworthPolynomial[order, cuttoff, variable, options];
 
order: filter order
cuttoff: butterworth cuttoff frequency

returns {poles,polynomial}

poles: list of stable poles of the Butterworth H(s)
polynomial: Buttworth polynomial in symbolic expression expressed in the variable supplied. The form of the polynomial depends on options used.

Possible values for options are 
form->\"linear\" returns H(s) in the form linear terms in the variable
form->\"quadratic\" returns H(s) in quadratic terms in the variable
If this option is not selected, then H(s) will be returned as one polynomial.

Example usage

Example 1
-----------

   Clear[s,form];
order=2;
cuttoff=.2 Pi;
{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s];
poles

{-0.44428829381583657`+0.44428829381583657` \[ImaginaryI],-0.44428829381583657`-0.44428829381583657` \[ImaginaryI]}

TraditionalForm@hs
\!\(\*FormBox[FractionBox[\(0.3947841760435742`\), \(\*SuperscriptBox[\(s\), \(2\)] + 0.8885765876316731`\\\ s + 0.39478417604357424`\)],
TraditionalForm]\)

Example 2
-----------

{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s,form\[Rule]\"linear\"];

poles
{-0.44428829381583657`+0.44428829381583657` \[ImaginaryI],-0.44428829381583657`-0.44428829381583657` \[ImaginaryI]}

TraditionalForm@hs
\!\(\*FormBox[\(\*FractionBox[\(0.4442882938158365`\\\ \[ImaginaryI]\), \(s + \((0.44428829381583657`\[InvisibleSpace] + 0.44428829381583657`\\\ \[ImaginaryI])\)\)] - \*FractionBox[\(0.4442882938158365`\\\ \[ImaginaryI]\), \(s + \((0.44428829381583657`\[InvisibleSpace] - 0.44428829381583657`\\\ \[ImaginaryI])\)\)]\),
TraditionalForm]\)

Example 3
-----------
{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s,form\[Rule]\"quadratic\"];
poles

{-0.44428829381583657`+0.44428829381583657` \[ImaginaryI],-0.44428829381583657`-0.44428829381583657` \[ImaginaryI]}

TraditionalForm@hs
\!\(\*FormBox[FractionBox[\(0.3947841760435742`\), \(\*SuperscriptBox[\(s\), \(2\)] + 0.8885765876316731`\\\ s + 0.39478417604357424`\)],
TraditionalForm]\)

"

(*----------------*)
(*                *)
(*----------------*)
butterd::"usage"="butterd[fs, variable, pass, stop, apass, astop, method, form];  

design a low pass IIR digital filter from Butterwoth analog filter.

fs   : sampling frequency in Hz
variable: the variable to use. example z
pass : passband corner frequency in Hz
stop : stopband corner frequency in Hz
apass : attenuation at passband in db, negative value
astop : attenuation at stopband in db, negative value
method: design method to use, can be  method->\"bilinear\", method->\"impulse\", method->\"matchedz\". If not 
given, then the default is \"bilinear\"
form: the form to return H(z) in. can be form->\"linear\" returns H(z) in the form 
linear terms in the variable, form->\"quadratic\" returns H(z) in quadratic 
terms in the variable. If this option is not selected, then H(z) will be returned as one polynomial.

Returns {order, cuttoff,digitalFilter, digitalFilterPoles, digitalFilterzeros, a, b}, 

order: Butterworth filter order found from the specifications
cuttoff: Butterworth cuttoff frequency found from the specifications
digitalFilter:  H(z)= B(z)/A(z), transfer function in the form selected above.
digitalFilterPoles:  locations of H(z) poles
digitalFilterzeros:  locations of H(z) zeros
a : the order+1 filter coefficients for A(z). Coefficients are listed in descending powers of z
b : the order+1 filter coefficients for B(z). Coefficients are listed in descending powers of z

EXAMPLE usage:

  Example1
------------

    test1
Clear[z];
fs=20000 (*hz*);
pass=2000 (* hz*);
stop=3000 (*hz*);
apass=-1 (*db*);
astop=-15 (*db*);
{filterOrder,\[CapitalOmega]c,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,form\[Rule]\"quadratic\"];
TraditionalForm@hz
\!\(\*FormBox[\(\((\(-\(0.8786598279577791`\\\ \*SuperscriptBox[\(z\), \(2\)]\)\) - 0.37458656177813776`\\\ z + 0.5040732661796413`)\)/\((1.`\\\ \*SuperscriptBox[\(z\), \(2\)] - 1.0105788810461338`\\\ z + 0.35827133770639796`)\) + \((0.05327834218253358`\\\ \*SuperscriptBox[\(z\), \(2\)] - 0.1260013381027728`\\\ z - 0.1792796802853064`)\)/\((1.`\\\ \*SuperscriptBox[\(z\), \(2\)] - 1.2686468043871897`\\\ z + 0.7051282432185233`)\) + \((0.8261193057058389`\\\ \*SuperscriptBox[\(z\), \(2\)] + 0.5806131391232774`\\\ z - 0.24550616658256155`)\)/\((1.`\\\ \*SuperscriptBox[\(z\), \(2\)] - 0.9043660641139339`\\\ z + 0.2155157075998965`)\)\),
TraditionalForm]\)
poles
{0.6343234021935948+0.550238189012782 I,
0.5052894405230669+0.320864331146801 I,
0.4521830320569668+0.10510096631175422 I,
0.4521830320569668-0.10510096631175422 I,
0.5052894405230669 -0.320864331146801 I,
0.6343234021935948-0.550238189012782 I}
zeros
{-1,-1,-1,-1,-1,-1}

  Example 2
------------
{filterOrder,\[CapitalOmega]c,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,method\[Rule]\"impulse\"];
TraditionalForm@hz
\!\(\*FormBox[\(\((12.619276514087685`\\\ \*SuperscriptBox[\(z\), \(5\)] + 202.07004065216458`\\\ \*SuperscriptBox[\(z\), \(4\)] + 322.8682701355137`\\\ \*SuperscriptBox[\(z\), \(3\)] + 82.01389599025786`\\\ \*SuperscriptBox[\(z\), \(2\)] + 2.0650372218951816`\\\ z)\)/\((\*SuperscriptBox[\(z\), \(6\)] - 3.3635196107885204`\\\ \*SuperscriptBox[\(z\), \(5\)] + 5.068420161781825`\\\ \*SuperscriptBox[\(z\), \(4\)] - 4.275864216158864`\\\ \*SuperscriptBox[\(z\), \(3\)] + 2.1066205743828084`\\\ \*SuperscriptBox[\(z\), \(2\)] - 0.5706492537421755`\\\ z + 0.06607428351012314`)\)\),
TraditionalForm]\)
poles
{0.6485799325392105`\[InvisibleSpace]+0.5236709777967435` \[ImaginaryI],0.5345537369865101`\[InvisibleSpace]+0.2901159614276238` \[ImaginaryI],0.49862613586853954`\[InvisibleSpace]+0.09176688744135254` \[ImaginaryI],0.49862613586853954`\[InvisibleSpace]-0.09176688744135254` \[ImaginaryI],0.5345537369865101`\[InvisibleSpace]-0.2901159614276238` \[ImaginaryI],0.6485799325392105`\[InvisibleSpace]-0.5236709777967435` \[ImaginaryI]}
zeros
{\[Infinity],\[Infinity],\[Infinity],\[Infinity],\[Infinity],\[Infinity]}

  Example 3
------------
{filterOrder,\[CapitalOmega]c,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,form\[Rule]\"linear\"];
TraditionalForm@hz
\!\(\*FormBox[\(\*FractionBox[\(3128.118552949461`\[InvisibleSpace] - 5418.060265807306`\\\ \[ImaginaryI]\), \(1 - \*FractionBox[\(0.6055598764724284`\[InvisibleSpace] - 0.5530645361185251`\\\ \[ImaginaryI]\), \(z\)]\)] - \*FractionBox[\(23348.59474341246`\), \(1 - \*FractionBox[\(0.4983854027120604`\[InvisibleSpace] + 0.2999718179941991`\\\ \[ImaginaryI]\), \(z\)]\)] - \*FractionBox[\(23348.59474341246`\), \(1 - \*FractionBox[\(0.4983854027120604`\[InvisibleSpace] - 0.2999718179941991`\\\ \[ImaginaryI]\), \(z\)]\)] + \*FractionBox[\(20220.476190462996`\[InvisibleSpace] - 35022.89211511871`\\\ \[ImaginaryI]\), \(1 - \*FractionBox[\(0.46770597726956914`\[InvisibleSpace] + 0.09398839450941299`\\\ \[ImaginaryI]\), \(z\)]\)] + \*FractionBox[\(20220.476190462996`\[InvisibleSpace] + 35022.89211511871`\\\ \[ImaginaryI]\), \(1 - \*FractionBox[\(0.46770597726956914`\[InvisibleSpace] - 0.09398839450941299`\\\ \[ImaginaryI]\), \(z\)]\)] + \*FractionBox[\(3128.118552949461`\[InvisibleSpace] + 5418.060265807306`\\\ \[ImaginaryI]\), \(1 - \*FractionBox[\(0.6055598764724284`\[InvisibleSpace] + 0.5530645361185251`\\\ \[ImaginaryI]\), \(z\)]\)]\),
TraditionalForm]\)
zeros
{-1,-1,-1,-1,-1,-1}

"

(*-----------------------------*)
(*                             *)
(*-----------------------------*)
getPositiveAndNegativeTerms::"usage"="{neg,pos}=getPositiveAndNegativeTerms[expr] 
Takes an expression expr and returns the terms in it which have negative sign and the terms which have postive sign

EXAMPLE

In[24]:= expr=3-c+b+d-h
{neg,pos}=getPositiveAndNegativeTerms[expr]

Out[25]= {{-c,-h},{3,b,d}}

Thanks to a post by Albert Retey on Math news group for the above idea of a solution to a question posted, which I saw on sept 9, 2010.
"

(*----------------*)
(*                *)
(*----------------*)
convertToLCD::"usage"="convertToLCD[mat]

convert matrix to least common denominator form.

Examples: 

convertToLCD[{{5/6, 1/4}, {1/7, 1/8}}] 
convertToLCD[{{5/6, 1/4}}] 

{{140/168,42/168},{24/168,21/168}}
{{10/12,3/12}}

Note that the result contains a HoldForm in it hence can't be used 
for computation as is without first releasing the hold"

(*----------------*)
(*                *)
(*----------------*)
checkForSingularity::"usage"="checkForSingularity[p, q, at, var]

Classifies singularity in ODE in the form y''+p(x) y' + q(x) y = 0

at is list of the points to check singularity at. The var is the
indepenent variable symbol, such as x. It returns a matrix, where
each row is of form {point, p limit, q limit}

If any one of the p limit or q limit is infinity, it means the singularity
is not removable. If the limits are finite, it means this is removable
singularity.

Examples: 

ClearAll[c, a, b, x];
p=(c - (a + b + 1) x)/(x (1 - x));
q=(-a b)/(x (1 - x));
at={0,1,Infinity};

m = checkForSingularity[p,q ,at, x];
Grid[Join[{{\"point\", \"limit x p(x)\", \"limit x^2 q(x)\"}}, m], 
 Frame -> All] 
 {{-1, 1, -(mu^2/4)}, {1, 1, -(mu^2/4)}, {\[Infinity], 0, \[Infinity]}}

Another example

ClearAll[c, a, b, x];
p=-2 x/(1 - x^2);
q=(lam + 40 (1 - x^2) - (mu^2)/(1 - x^2))/(1 - x^2);
at={0,1,Infinity};

m = checkForSingularity[p,q ,at, x];

{{-1, 1, -(mu^2/4)}, {1, 1, -(mu^2/4)}, {\[Infinity], 0, \[Infinity]}}

"


(*==========================================================*)
(*==========================================================*)
(*==========================================================*)
(*==========================================================*)
Begin["`Private`"]

getStateGainVector[a_, b_, eigs_List, debug_] :=
    Module[ {nA, cA, nB, cB, A0, B0, sysOriginal,
      cOriginal, sysTarget, cTarget, transformationT, gains, k, aClosed, pClosed, s, pClosedCoefficnets,
      pDesign, pDesignCoefficnets, gainValues},
      messageHandler = If[ Last[#], Abort[]]&;
      Internal`AddHandler["Message", messageHandler];
      General::error = "`1`";

      {nA, cA} = Dimensions[a];
      If[ nA != cA, Message[getStateGainVector::error, "Matrix A must be square"]; Return[] ];

      {nB, cB} = Dimensions[b];
      If[ nB != nA, Message[getStateGainVector::error, "Matrix B must have same rows as A"]; Return[] ];

      If[ cB != 1,
        Message[getStateGainVector::error, "B must be one column matrix. Only single input is supported"]; Return[]];

      If[ Length[eigs] != nA,
        Message[getStateGainVector::error, "number of design eigenvalues much be the same as size of A"]; Return[] ];
      If[ debug,
        Print["A=", MatrixForm[a]]
      ];
      If[ debug,
        Print["B=", MatrixForm[b]]
      ];

      (*check if system is controllable*)
      sysOriginal = StateSpaceModel[{a, b}];
      If[ Not[ControllableModelQ[sysOriginal]],
        Message[getStateGainVector::error, "System must be controllable"];
        Return[]
      ];

      (*convert to controllable form (we call this the target system)*)
      {A0, B0} = nma`getControllableForm[a, b];
      If[ debug,
        Print["A in controllable format=", MatrixForm[A0]]
      ];
      If[ debug,
        Print["B in controllable format=", MatrixForm[B0]]
      ];

      (*Obtain controllability matrix for original system*)
      cOriginal = ControllabilityMatrix[sysOriginal];
      If[ debug,
        Print["Controllability matrix for original system=", MatrixForm[cOriginal]]
      ];

      (*Obtain controllability matrix for target system*)
      sysTarget = StateSpaceModel[{A0, B0}];
      cTarget = ControllabilityMatrix[sysTarget];
      If[ debug,
        Print["Controllability matrix for target system=", MatrixForm[cTarget]]
      ];

      (*Obtain T,the transformation matrix*)
      transformationT = cTarget.Inverse[cOriginal];
      If[ debug,
        Print["transformation T matrix = ", MatrixForm[transformationT]]
      ];

      (*define gains*)
      gains = Array[k, Length[eigs]];

      (*generate the closed loop state feedback A matrix using the traget A,B system*)
      aClosed = A0 + B0.{gains};
      If[ debug,
        Print["Closed loop state feedback = ", MatrixForm[nma`str@aClosed]]
      ];

      (*Find the CharacteristicPolynomial of the above matrix*)
      pClosed = CharacteristicPolynomial[aClosed, s];
      If[ debug,
        Print["Closed loop Characteristic Polynomial = ", nma`str[Simplify@pClosed]]
      ];

      (*Extract the coefficients so we can compare them to the design polynomial below*)
      pClosedCoefficnets = CoefficientList[pClosed, s];
      If[ debug,
        Print["Closed loop Characteristic Polynomial coefficients = ", nma`str@pClosedCoefficnets]
      ];

      (*Generate the design polynomial from the desired pole locations*)
      pDesign = Expand[Times @@ ((s - #)& /@ eigs)];
      If[ debug,
        Print["Design Characteristic Polynomial = ", nma`str[Simplify@pDesign]]
      ];

      (*Extract the coefficients of the design polynomial*)
      pDesignCoefficnets = CoefficientList[pDesign, s];
      If[ debug,
        Print["design Characteristic Polynomial coefficients = ", pDesignCoefficnets]
      ];

      (*Solve for the gains by comparing coefficients of design polynomial with the closed loop polynomial*)

      (*make sure a0 coefficient have same sign*)
      If[ Sign[Last@pDesignCoefficnets] < 0,
        If[ Sign[Last@pClosedCoefficnets] > 0,
          pDesignCoefficnets = -pDesignCoefficnets
        ],
        If[ Sign[Last@pClosedCoefficnets] < 0,
          pClosedCoefficnets = -pClosedCoefficnets
        ]
      ];
      gainValues = gains /. First@Solve[pDesignCoefficnets == pClosedCoefficnets, gains];
      If[ debug,
        Print["gains in target space are = ", gainValues]
      ];

      (*Convert the gain vector to original space using the transformation T found above*)
      gainValues = gainValues.transformationT;
      If[ debug,
        Print["gains in original space are = ", gainValues]
      ];
      gainValues
    ];



getControllableForm[a_, b_] :=
    Module[ {sys, tf},
      messageHandler = If[ Last[#],
        Abort[]
      ]&;
      Internal`AddHandler["Message", messageHandler];
      sys = StateSpaceModel[{a, b}];
      tf = TransferFunctionModel[sys];
      sys = StateSpaceModel[tf, StateSpaceRealization -> "Controllable"];
      {Normal[sys][[1]], Normal[sys][[2]]}
    ];


generateLTVObservabilityMatrix[a_List, c_List, t_] :=
    Module[ {nA, mA, nC, mC, matrix, mLast, k},
      General::error = "`1`";
      {nA, mA} = Dimensions[a];
      {nC, mC} = Dimensions[c];
      If[ nA != mC,
        Message[generateLTVObservabilityMatrix::"error", "number of rows in A not same as number of rows in B"];
        Return[]
      ];
      If[ nA != mA,
        Message[generateLTVObservabilityMatrix::"error", "A must be square matrix"];
        Return[]
      ];
      mLast = c;
      matrix = mLast;
      Do[
        mLast = mLast.a + D[mLast, t];
        matrix = Join[matrix, mLast]
        ,
        {k, 0, nA - 2}
      ];
      matrix
    ];



generateLTVControllabilityMatrix[a_List, b_List, t_] :=
    Module[ {nA, mA, nB, mB, matrix, mLast, k},
      General::error = "`1`";
      {nA, mA} = Dimensions[a];
      {nB, mB} = Dimensions[b];
      If[ nA != nB,
        Message[generateLTVControllabilityMatrix::"error",
          "number of rows in A not same as number of rows in B"];
        Return[]
      ];
      If[ nA != mA,
        Message[generateLTVControllabilityMatrix::"error", "A must be square matrix"];
        Return[]
      ];
      mLast = b;
      matrix = mLast;
      Do[
        mLast = -a.mLast + D[mLast, t];
        matrix = Transpose@Join[Transpose[matrix], Transpose@mLast]
        ,
        {k, 0, nA - 2}
      ];
      matrix
    ];


displayHurwitz[p_List, s_] :=
    Module[ {z, i, roots, n, leading},
      plot[data_] :=
          (
            ListPlot[{Re[#], Im[#]}& /@ data, AxesOrigin -> {0, 0}, ImagePadding -> 40, AspectRatio -> 1, Frame -> True,
              FrameLabel -> {{Im, None}, {Re, "complex plane"}}, PlotStyle -> Directive[Red, PointSize[.06]]]
          );
      n = Length[CoefficientList[First@p, s]] - 1;
      leadingLabel = Table[Subscript["\[CapitalDelta]", i], {i, 1, n}];
      data = {
        z = hurwitz[#, s];
        Grid[{{
          Style[TraditionalForm@#, 12]
        },
          {MatrixForm[z]}
        }, Frame -> All, FrameStyle -> Gray, Spacings -> {.1, .5}],
        leading = leadingDet[z];
        Grid[Transpose[{leadingLabel, leading}], Frame -> All],
        roots = s /. NSolve[# == 0, s];
        plot[N[roots]],
        Column[Re[roots]]}& /@ p;
      header = {"Hurwitz Matrix", "\!\(\*SubscriptBox[\(\[CapitalDelta]\), \(i\)]\)", "root locations",
        "Real part of roots"};
      data = Insert[data, header, 1];
      Grid[data, Frame -> All]
    ];


hurwitz[p_, s_] :=
    Module[ {n, c, t, e, o, cr, cc, debug = False},
      c = CoefficientList[p, s];
      If[ debug,
        Print["c=", c]
      ];
      n = Length[c];
      If[ debug,
        Print["n=", n]
      ];
      If[ c[[1]] < 0,
        c = -c
      ];
      e = If[ EvenQ[#],
        #,
        Sequence @@ {}
      ]& /@ Range[n];
      o = If[ OddQ[#],
        #,
        Sequence @@ {}
      ]& /@ Range[n];
      t = Table[0, {n - 1}, {n - 1}];
      cc = 1;
      cr = 1;
      done = False;
      While[Not[done],
        t[[cr, cc ;; cc + Length[e] - 1]] = c[[e]];
        If[ debug,
          Print["t=", t]
        ];
        cr++;
        If[ cr <= (n - 1),
          t[[cr, cc ;; cc + Length[o] - 1]] = c[[o]];
          cc++
        ];
        cr++;
        If[ cr > (n - 1),
          done = True
        ]
      ];
      t
    ];


leadingDet[h_] := Det[h[[1 ;; #, 1 ;; #]]]& /@ Range[Length[h]];


kharitonovPoly[w_List, s_] :=
    Module[ {n = Length[w], p = Table[0, {4}], z, i},
      
      z = 2;
      data = Flatten[Table[ If[ z == 2,
        z = 1; {1, 1},
        z = 2; {2, 2}
      ], {i, 1, Ceiling[n / 2]}]];

      data = data[[1 ;; n]];
      p[[1]] = Total@Table[ s^(i - 1)w[[i, data[[i]]]], {i, 1, Length[data]}];

      z = 1;
      data = Flatten[Table[ If[ z == 2,
        z = 1; {1, 1},
        z = 2; {2, 2}
      ], {i, 1, Ceiling[n / 2]}]];

      data = data[[1 ;; n]];
      p[[2]] = Total@Table[ s^(i - 1)w[[i, data[[i]]]], {i, 1, Length[data]}];

      z = 1;
      data = Flatten[Table[ If[ z == 2,
        z = 1; {1, 2},
        z = 2; {2, 1}
      ], {i, 1, Ceiling[n / 2]}]];

      data = data[[1 ;; n]];
      p[[3]] = Total@Table[ s^(i - 1)w[[i, data[[i]]]], {i, 1, Length[data]}];

      z = 2;
      data = Flatten[Table[ If[ z == 2,
        z = 1; {1, 2},
        z = 2; {2, 1}
      ], {i, 1, Ceiling[n / 2]}]];

      data = data[[1 ;; n]];
      p[[4]] = Total@Table[ s^(i - 1)w[[i, data[[i]]]], {i, 1, Length[data]}];
      p
    ];


(*This function below is not mine, but obtained from *)
(*http://library.wolfram.com/infocenter/MathSource/5026/ *)
(*the Author is Peter Valk\[CapitalATilde]\.b3, I put it here so I can use it easily,*)
(*it does numerical inversion of Laplace transform*)
(*example calling *)
(*fun[s_]=(1/s) Exp[-Sqrt[ (s^2+ (37/100) s+1)/(s^2+ s+Pi)] ];*)
(*t0=100;*)
(*ft=FT[fun,t0];*)
(*------------------------------------------------------------*)

FT[F_, t_, M_ : 32] := Module[{np, r, S, theta, sigma}, np = Max[M, $MachinePrecision];
r = SetPrecision[2M / (5t), np];
S = r theta (Cot[theta] + I);
sigma = theta + (theta Cot[theta] - 1)Cot[theta];
(r / M)Plus @@ Append[Table[Re[Exp[t S](1 + I sigma)F[S]], {theta, Pi / M, (M - 1)Pi / M, Pi / M}], (1 / 2) Exp[r t] F[r]]];


matrix2DredBlackPositions[
  nRow_?(IntegerQ[#] && Positive[#] &),
  nCol_?(IntegerQ[#] && Positive[#] &),
  OptionsPattern[]
] :=
    Module[ {red, black, isInternal = OptionValue[internal], nrow, ncol, j},
      If[ OddQ[nCol],
        {
          If[ isInternal,
            {nrow = nRow - 2;
            ncol = nCol - 2},
            {nrow = nRow;
            ncol = nCol}
          ];
          red = Position[Partition[Table[(-1)^j, {j, Range[nrow * ncol]}], ncol], -1];
          black = Position[Partition[Table[(-1)^j, {j, Range[nrow * ncol]}], ncol], 1];
          If[ isInternal,
            {red += 1;
            black += 1}
          ];
        },
        {
          Throw["matrix2DredBlackPositions:: supports odd columns only at this time"]
        }
      ];
      {red, black}
    ];
Options[matrix2DredBlackPositions] = {internal -> True};


matrix2DLexOrdering[
  mat_?(MatrixQ[#]&)
] :=
    Module[ {nRow, nCol, s, A, B},
      If[ Length[Dimensions[mat]] != 2,
        Throw["matrix2DLexOrdering::  matrix must be 2D"]
      ];
      {nRow, nCol} = Dimensions[mat];
      s = Outer[List, Range[nRow], Range[nCol]];
      A = Reverse[Reverse[s, 3]];
      B = Reverse@Partition[Table[i, {i, nRow * nCol}], nCol];
      {A, B}
    ];


centerImageData[
  im_?(MatrixQ[#]&)
] :=
    Module[ {nRow, nCol, indx, d, i, j},
      If[ Length[Dimensions[im]] != 2,
        Throw["centerImageData:: image matrix must be 2D"]
      ];
      {nRow, nCol} = Dimensions[im];
      d = Table[0, {i, nRow}, {j, nCol}];
      indx = Flatten[Outer[List, Range[1, nRow], Range[1, nCol]], 1];
      Map[(d[[#[[1]], #[[2]]]] = (-1)^(#[[1]] + #[[2]]) * im[[#[[1]], #[[2]]]])&, indx];
      d
    ];


diagonalBlockMatrix[
  block_?(MatrixQ[#]&),
  nBlocks_?(IntegerQ[#] && Positive[#] &)
] :=
    Module[ {nRow, nCol},
      {nRow, nCol} = Dimensions[block];
      If[ nRow != nCol,
        Throw["diagonalBlockMatrix:: Matrix must be square"]
      ];
      SparseArray[Band[{1, 1}] -> ConstantArray[block, {nBlocks}]]
    ];


poissonMatrixOnSquareDirichlet[n_?(IntegerQ[#] && Positive[#] &),
  m_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {block, sizeOfBlock, nBlocks, a, i, sub, r},
      If[ m != 5,
        Throw["poissonMatrixOnSquareDirichlet:: m = 5 is only value currently supported"]
      ];

      (* first make the block*)
      r = Table[-4, {i, n}];
      block = DiagonalMatrix[r, 0];
      r = Table[1, {i, n - 1}];
      block += DiagonalMatrix[r, 1];
      block += DiagonalMatrix[r, -1];
      sizeOfBlock = Length[block];
      nBlocks = n;
      a = Normal@diagonalBlockMatrix[block, nBlocks];
      sub = Table[1, {i, (n - 1)n}];
      a = a + DiagonalMatrix[sub, n];
      a = a + DiagonalMatrix[sub, -n]
    ];


diagonalDominantMatrixQ[mat_?(MatrixQ[#]&), OptionsPattern[]] :=
    Module[ {isStrict = OptionValue[strict], nRow, nCol, d, i, j, s, result = False, keepChecking},
      If[ isStrict =!= True && isStrict =!= False,
        Throw["diagonalDominantMatrixQ:: invalid option detected"]
      ];
      {nRow, nCol} = Dimensions[mat];
      If[ nRow != nCol,
        Throw["diagonalDominantMatrixQ:: Matrix must be square"]
      ];
      keepChecking = True;
      For[i = 1, (i <= nRow && keepChecking), i++,
        {
          d = Abs[mat[[i, i]]];
          s = Sum[Abs[mat[[i, j]]], {j, nCol}];
          If[ isStrict,
            If[ d > (s - d),
              {result = True;
              keepChecking = False}
            ],
            If[ d >= (s - d),
              {result = True;
              keepChecking = False}
            ]
          ];

        }];
      result
    ];
Options[diagonalDominantMatrixQ] = {strict -> False};


numItse[v_?(NumberQ[#]&), (*the number to format*)
  s1_?(IntegerQ[#] && Positive[#] &), (* total number of digits to use, including the decimal digits*)
  s2_?(IntegerQ[#] && Positive[#] &), (*number of digits to display to the right of the decimal point*)
  fontSize_?(IntegerQ[#] && Positive[#] &) (*font size to use for string formating*)
] :=
    Module[ {},
      Style[ScientificForm[N[v, $MachinePrecision], {s1, s2},
        NumberFormat -> (Row[{#1, "e", #3}]&), NumberPadding -> {"", "0"}], fontSize]
    ];


linspace[x1_?(NumericQ[#]&),
  x2_?(NumericQ[#]&),
  n_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {h},
      If[ n < 2,
        Return[x2]
      ];
      h = (x2 - x1) / (n - 1);
      Range[x1, x2, h]
    ];



str[expr_] := Module[{},
  StringReplace[ToString[expr, FormatType -> TraditionalForm],
    c : LetterCharacter ~~ "$" ~~ DigitCharacter.. :> c]
];


quantizationLevel[x_?(NumericQ[#]&),
  mp_?(NumericQ[#] && Positive[#] &),
  nBits_?(IntegerQ[#] && Positive[#] &)] := Module[{nStep, stepSize},

  stepSize = (2 mp) / 2^nBits;
  nStep = If[x < 0,
    Round[(mp - Abs[x]) / stepSize],
    Round[x / stepSize]
  ]
];


toTwoComplement[x_?(NumericQ[#]&), (*analog value*)
  mp_?(NumericQ[#] && Positive[#] &), (*maximum allowed level*)
  nBits_?(IntegerQ[#] && Positive[#] &) (*number of bits available for coding*)
] :=
    Module[ {q, delta, b = nBits - 1, code},
      delta = mp / 2^b;  (*minumum difference between each level*)
      q = Round[Abs[x] / delta]; (*quantization level*)

      (*now determine the code of x based on two complement logic *)
      If[ x >= -delta / 2,
        {
          If[ q >= (2^b) - 1,
            q = (2^b) - 1
          ];
          code = q;
        },
        {
          If[ q >= (2^b) - 1,
            q = (2^b)
          ];
          code = 2^nBits - q;
        }
      ];
      {N@delta, N[Sign[x] * q * delta], code}
    ];


toOffsetBinary[x_?(NumericQ[#]&),
  mp_?(NumericQ[#] && Positive[#] &),
  nBits_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {q, delta, b = nBits - 1, code},
      delta = mp / 2^b;
      q = Round[Abs[x] / delta];
      If[ x >= -delta / 2,
        {
          If[ q >= (2^b) - 1,
            q = (2^b) - 1
          ];
          code = 2^b + q;
        },
        {
          If[ q >= (2^b) - 1,
            q = (2^b)
          ];
          code = 2^b - q;
        }
      ];
      {N@delta, N[Sign[x] * q * delta], code}
    ];


toOnesComplement[x_?(NumericQ[#]&),
  mp_?(NumericQ[#] && Positive[#] &),
  nBits_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {code, delta, b = nBits - 1, q},
      delta = mp / 2^b;
      q = Round[Abs[x] / delta];
      If[ q >= (2^b) - 1,
        q = (2^b) - 1
      ];
      If[ x > 0,
        code = q
        ,
        code = (2^nBits - 1) - q
      ];
      {N@delta, N[Sign[x] * q * delta], code}
    ];


toSignMagnitude[x_?(NumericQ[#]&),
  mp_?(NumericQ[#] && Positive[#] &),
  nBits_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {code, delta, b = nBits - 1, q},
      delta = mp / 2^b;
      q = Round[Abs[x] / delta];
      If[ q >= (2^b) - 1,
        q = (2^b) - 1
      ];
      If[ x >= 0,
        code = q
        ,
        code = q + 2^b
      ];
      {N@delta, N[Sign[x] * q * delta], code}
    ];


formatBinary[num_?(IntegerQ[#] && (# >= 0) &), nBits_?(IntegerQ[#] && Positive[#] &)] :=
    Module[ {s = "", r, i, j, a, nRow, nCol},
      r = IntegerDigits[num, 2];
      r = PadLeft[r, nBits];
      a = Partition[r, If[ nBits < 4,
        nBits,
        4
      ]];
      {nRow, nCol} = Dimensions[a];
      Do[
        Do[

          If[ i > 1 && j == 1,
            s = StringJoin[s, " " ];
          ];
          s = StringJoin[s, ToString[ a[[i, j]]  ]];

          , {j, 1, nCol}
        ],
        {i, 1, nRow}
      ];
      s
    ];


getPhaseOfPulseTrainFourierCoeff[
  delay_?NumericQ,
  period_?NumericQ,
  range_?NumericQ,
  dutyCycle_?NumericQ,
  nCoeff_Integer] :=
    Module[ {data, n, t0},
      t0 = Mod[delay, period];
      data = {{0, 0}};
      data = Append[List@Reverse@Table[{n,
        getPhase[n, Sinc[Pi n dutyCycle], t0, period]}, {n, -1, -(nCoeff - 1), -1}], data];
      data = Append[data, Table[{n,
        getPhase[n, Sinc[Pi n dutyCycle], t0, period]}, {n, 1, nCoeff - 1, 1}]];
      Chop@Flatten[data, 1]
    ];


getMagnitudeOfPulseTrainFourierCoeff[
  t0_?NumericQ,
  period_?NumericQ,
  range_?NumericQ,
  dutyCycle_?NumericQ,
  nCoeff_Integer,
  h_?(NumberQ[#] && Positive[#]&)] :=
    Module[ {data, n},
      data = {{0, h * dutyCycle}};

      data = Append[List@Reverse@Table[{n,
        Abs[h * dutyCycle * Sinc[Pi * n * dutyCycle] * Exp[-I * 2 * Pi / period * n * t0]]},
        {n, -1, -(nCoeff - 1), -1}], data];

      data = Append[data,
        Table[{n, Abs[h * dutyCycle * Sinc[Pi * n * dutyCycle] * Exp[-I * 2 * Pi / period * n * t0]]},
          {n, 1, nCoeff - 1, 1}]];

      Chop@Flatten[data, 1]
    ];


makePulseTrain[delay_?(NumberQ[#] && # >= 0&), (*delay of a pulse from time=0*)
  period_?(NumberQ[#] && Positive[#]&), (* the time between each pulse*)
  range_?(NumberQ[#] && Positive[#]&), (* the time duration to use*)
  tao_?(NumberQ[#] && Positive[#]&), (* width of each pulse*)
  h_?(NumberQ[#] && Positive[#]&)        (* height of the pulse*)
] := Module[{allPulses},
  getForwardPulses[] :=
      Module[ {n, z, pulse, forwardPulses = {}, more, t0},
        t0 = Mod[delay, period];
        n = 0;
        more = True;
        While[more,
          {
            z = t0 + n * period;
            If[ z > range,
              more = False,
              {
                pulse = Line[{ {z - tao / 2, 0}, {z - tao / 2, h}, {z + tao / 2, h}, {z + tao / 2, 0}}];
                forwardPulses = Append[forwardPulses, pulse];
                n = n + 1;
              }
            ]
          }
        ];
        forwardPulses
      ];

  getBackPulses[] :=
      Module[ {n, z, pulse, backPulses = {}, more, t0},
        t0 = Mod[delay, period];
        n = 0;
        more = True;
        While[more,
          {
            z = -(period - t0) - n * period;
            If[ Abs[z] > range,
              more = False,
              {
                pulse = Line[{ {z + tao / 2, 0}, {z + tao / 2, h}, {z - tao / 2, h}, {z - tao / 2, 0}}];
                backPulses = Append[backPulses, pulse];
                n = n + 1;
              }
            ]
          }
        ];
        backPulses
      ];

  allPulses = getForwardPulses[];
  allPulses = Append[allPulses, getBackPulses[]]
];


fcPulseTrain[
  n_?(IntegerQ[#]&), (* nth coeff *)
  h_?(NumberQ[#] && Positive[#]&), (*hight of pulse *)
  width_?(NumberQ[#] && Positive[#]&), (*width of pulse *)
  period_?(NumberQ[#] && Positive[#]&) (*period of pulse *)
] := Module[ {dutyCycle = width / period},
  h * dutyCycle * Sinc[Pi * n * dutyCycle]
];


numIt[v_?(NumberQ[#]&), (*the number to format*)
  s1_?(IntegerQ[#] && Positive[#] &), (* total number of digits to use, including the decimal digits*)
  s2_?(IntegerQ[#] && Positive[#] &), (*number of digits to display to the right of the decimal point*)
  fontSize_?(IntegerQ[#] && Positive[#] &) (*font size to use for string formating*)
] := Module[ {},
  Style[ToString[AccountingForm[Chop[v], {s1, s2},
    NumberPadding -> {" ", "0"}, NumberSigns -> {"-", ""}]], fontSize]
];


butterToLowPass[h_,
  \[CapitalOmega]c_?(NumericQ[#] && Positive[#] &),
  s_] := Module[{},
  h /. s -> (1 / \[CapitalOmega]c s)
];


butterToHighPass[h_, \[CapitalOmega]c_?(NumericQ[#] && Positive[#] &), s_] := Module[{},
  h /. s -> \[CapitalOmega]c / s
];


butterToBandPass[h_,
  omega1_?(NumericQ[#] && Positive[#] &),
  omega2_?(NumericQ[#] && Positive[#] &),
  s_] := Module[{},
  If[omega2 < omega1,
    Throw["butterToBandPass: check on \[CapitalOmega]2<\[CapitalOmega]1"]
  ];
  h /. s ->( (s^2 + omega1 omega2)/(s(omega2 - omega1)) )
];


butterToBandStop[h_,
  omega1_?(NumericQ[#] && Positive[#] &),
  omega2_?(NumericQ[#] && Positive[#] &),
  s_] := Module[{},
  If[omega2 < omega1,
    Throw["butterToBandStop: check on \[CapitalOmega]2<\[CapitalOmega]1"]
  ];
  h /. s -> (s(omega2 - omega1)) / (s^2 + omega1 omega2)
];


getButterworthPolynomial[
  filterOrder_?(IntegerQ[#] && Positive[#] &), (*filter order, N *)
  omegac_?(NumericQ[#] && Positive[#]&), (*cutoff freq. must be less than Pi*)
  s_ , (*variable to use *)
  opts : OptionsPattern[]
  ] := Module[ {poles, gain, coeff, hs, form = OptionValue[form]},  
      If[omegac >= Pi,
        Throw["getButterworthPolynomial: cut off frequency can not be larger than Pi"]
      ];
      poles = getButterworthStablePoles[filterOrder, omegac];
      gain  = Chop@N@getButterworthGain[poles];
      coeff = getPartialFractionsCoeff[poles, gain];
      If[ Order[form, "linear"] == 0,
        hs = Total@Table[ coeff[[n]] / (s - poles[[n]]), {n, 1, Length[poles]}],
        hs = getButterworthQuadraticTerms[poles, coeff, s];
        If[ Order[form, "quadratic"] == 0,
          hs = Total@hs,
          hs = FullSimplify@Total@hs;
          hs = Chop@Expand[Numerator@hs] / Chop@Expand@Denominator@hs
        ]
      ];
      {poles, hs}
    ];

Options[getButterworthPolynomial] = {form -> "default"};


butterd[
  fs_?(NumberQ[#] && Positive[#]&), (* sampling frequency in Hz*)
  var_, (*symbol for H(z) *)
  pass_?(NumberQ[#] && Positive[#]&), (* passband corner frequency in Hz*)
  stop_?(NumberQ[#] && Positive[#]&), (* stopband corner frequency in Hz*)
  deltap_?(NumberQ[#] && Negative[#]&), (* attenuation at passband in db*)
  deltas_?(NumberQ[#] && Negative[#]&), (* attenuation at stopband in db*)
  OptionsPattern[]
] :=
    Module[ {T = 1 / fs, As, Ap, filterOrder, omegac, omegapass, omegastop, wpass, wstop,
      method = OptionValue[method], form = OptionValue[form], hz, hspoles, coeff, gain, hs, s, myhz,
      maxOrder, DEBUG = False, c, hzpoles, hzzeros, numeratorCoeff, denominatorCoeff},
      If[ DEBUG,
        Print["enter butterd"]
      ];
      wpass = (pass / fs) 2 Pi;
      wstop = (stop / fs) 2 Pi;
      If[ Order[method, "bilinear"] == 0,
        omegapass = N[2 / T Tan[wpass / 2]];
        omegastop = N[2 / T Tan[wstop / 2]],
        omegapass = N[ wpass / T];
        omegastop = N[ wstop / T]
      ];
      As = N[10^(-deltas / 10)];
      Ap = N[10^(-deltap / 10)];
      filterOrder = 1 / 2 ((Log[10, Ap - 1] - Log[10, As - 1]) / (Log[10, omegapass] - Log[10, omegastop]));
      filterOrder = Ceiling[filterOrder];
      If[ DEBUG,
        Print["filter order = ", filterOrder]
      ];
      omegac = N[If[ Order[method, "bilinear"] == 0,
        omegastop / 10^(1 / (2  filterOrder) Log[10, As - 1]),
        omegapass / 10^(1 / (2  filterOrder) Log[10, Ap - 1])
      ]];
      If[ DEBUG,
        Print["omegac = ", omegac]
      ];
      hspoles = getButterworthStablePoles[filterOrder, omegac];
      If[ DEBUG,
        Print["hspoles = ", hspoles]
      ];
      gain = Chop@N@getButterworthGain[hspoles];
      If[ DEBUG,
        Print["gain = ", gain]
      ];
      coeff = getPartialFractionsCoeff[hspoles, gain];
      If[ DEBUG,
        Print["coeff = ", coeff]
      ];
      If[ Order[method, "impulse"] == 0,
        If[ DEBUG,
          Print["in impulse method"]
        ];
        hz = getHzImpulseInvariance[hspoles, coeff, fs, var];,
        If[ DEBUG,
          Print["in bilinear method"]
        ];
        hs = getHsFromPolesSecondOrder[hspoles, s, coeff];
        If[ DEBUG,
          Print["hs=", hs]
        ];
        hz = getHzBilinear[hs, fs, s, var];
        If[ DEBUG,
          Print["hz=", hz]
        ]
      ];
      myhz = FullSimplify@Total@hz;
      myhz = Chop@Expand[Numerator@myhz] / Chop@Expand@Denominator@myhz;
      numeratorCoeff = CoefficientList[Numerator[myhz], var];
      denominatorCoeff = CoefficientList[Denominator[myhz], var];
      If[ Order[form, "linear"] == 0,
        If[ DEBUG,
          Print["in linear form"]
        ];
        hz = Table[ (coeff[[i]]) / (1 - var^-1 * Exp[hspoles[[i]] T]) , {i, 1, Length[hspoles]}];
        hz = Total@hz,
        If[ Order[form, "quadratic"] == 0,
          If[ DEBUG,
            Print["in quadratic form"]
          ];
          myhz = Table[0, {i, Length[hz]}];
          Do[
            myhz[[i]] = Expand@FullSimplify[Numerator[hz[[i]]]] /
                Expand@FullSimplify[Denominator[hz[[i]]]];
            If[ DEBUG,
              Print["myhz[[", i, "]=", myhz[[i]]]
            ];
            maxOrder = Exponent[ Denominator[hz[[i]]] , var];
            If[ DEBUG,
              Print["myhz[[", i, "], max order of denominator is ", maxOrder]
            ];
            maxOrder = Exponent[ Numerator[hz[[i]]] , var];
            If[ DEBUG,
              Print["myhz[[", i, "], max order of Numerator is ", maxOrder]
            ];
            c = CoefficientList[Denominator[hz[[i]]], var];
            If[ DEBUG,
              Print["c=", c]
            ];
            If[ Length[c] == 3,
              c = c[[3]];
              If[ DEBUG,
                Print["c[3]=", c]
              ];
              myhz[[i]] = Expand@FullSimplify[Numerator[myhz[[i]]] / c] /
                  Expand@FullSimplify[Denominator[myhz[[i]]] / c];
              If[ DEBUG,
                Print["myhz[", i, "] now is=", myhz[[i]] ]
              ];
            ]
            , {i, 1, Length[hz]}
          ];
          hz = Total@myhz;,
          hz = myhz
        ]
      ];
      hzpoles = Chop@getHzPoles[hspoles, method, fs];
      If[ Order[method, "impulse"] == 0,
        If[ DEBUG,
          Print["method is impulse"]
        ];
        hzzeros = Table[Infinity, {i, 1, Length[hzpoles]}];
        If[ DEBUG,
          Print["hzzeros=", hzzeros]
        ];,
        If[ Order[method, "bilinear"] == 0,
          hzzeros = Table[-1, {i, 1, Length[hzpoles]}],
          hzzeros = Table[-1, {i, 1, Length[hzpoles]}];
        ]
      ];
      {filterOrder, omegac, hz, hzpoles, hzzeros, numeratorCoeff, denominatorCoeff}
    ];
Options[butterd] = {method -> "bilinear", form -> "default"};


(*--------------------------------------------*)
(*  Need it to handle DirectDelta derivative  *)
(*--------------------------------------------*)
plotFourierTransform[ft_, w_, from_, to_, ymin_, ymax_, arrowhead_] :=
    Module[ {any, any1, terms, e, amplitude, shift, impulsesFound = 0, DEBUG = False},
      terms = Table[ft[[i]], {i, 1, Length[ft]}];
      Do[
        e = terms[[i]];
        If[ DEBUG,
          Print["e=", e]
        ];
        shift = e /. (any1_. DiracDelta[any_]) -> any;
        If[ UnsameQ[shift, e], (*Impulse found *)
          {
            If[ DEBUG,
              Print["found impulse, shift=", shift]
            ];
            impulsesFound = impulsesFound + 1;
            amplitude = e /. (a_. DiracDelta[any_]) -> a;
            If[ DEBUG,
              Print["amplitude=", amplitude]
            ];
            shift = shift /. (any_. + w ) -> any;
            If[ DEBUG,
              Print["shift now is =", shift]
            ];
            amplitude = amplitude /. w -> (-shift);
            If[ DEBUG,
              Print["amplitude now after subs is=", amplitude]
            ];
            If[ SameQ[Head[amplitude], Complex],
              {
                If[ DEBUG,
                  Print["found complex amplitude="]
                ];
                terms[[i]] = Graphics[{
                  {Red, Arrowheads[arrowhead], Arrow[{  {-shift, 0}, {-shift, Re[amplitude]}  }]},
                  {Red, Arrowheads[arrowhead], Arrow[{  {-shift, 0}, {-shift, Im[amplitude]}  }]}
                }];
                If[ DEBUG,
                  Print["term[", i, "]=", terms[[i]]]
                ];
              },
              {
                If[ DEBUG,
                  Print["not complex amplitude="]
                ];
                terms[[i]] = Graphics[
                  {Arrowheads[arrowhead], Arrow[{{-shift, 0}, {-shift, amplitude}}]}
                ];
                If[ DEBUG,
                  Print["term[", i, "]=", terms[[i]]]
                ];
              }
            ];
          },
          {
            terms[[i]] = Plot[terms[[i]], {w, from, to}];
          }
        ],
        {i, 1, Length[terms]}
      ];
      If[ DEBUG,
        Print["terms=", terms]
      ];
      If[ impulsesFound == Length[terms],
        {
          Show[Plot[0, {w, from, to}], terms, PlotRange -> {All, {ymin, ymax}}]
        },
        {
          Show[terms, PlotRange -> {All, {ymin, ymax}}]
        }
      ]
    (*If[impulsesFound==Length[terms],
        {
        {Plot[0,{w,from,to}],terms}
        }
        ,
        {
        terms
        }
        ]*)

    ];


getPositiveAndNegativeTerms[expr_] := Module[{neg, pos},
  neg = Cases[expr, _?Negative | Times[_?Negative, ___]];
  pos = Cases[expr, _?Positive | _Symbol | Times[_?Positive, ___]];
  {neg, pos}
];


(***************** Internal ************************)
(***************** Internal ************************)
(***************** Internal ************************)
(***************** Internal ************************)


getHzPoles[hsPoles_, method_, fs_] :=
    Module[ {T = 1 / fs, i},
      If[ Order[method, "bilinear"] == 0,
        Table[(1 + (T / 2)hsPoles[[i]]) / (1 - (T / 2)hsPoles[[i]]), {i, 1, Length[hsPoles]}],
        Table[Exp[T * hsPoles[[i]]], {i, 1, Length[hsPoles]}]
      ]
    ];


getHzBilinear[hsSecondOrderForm_List,
  sampling_,
  s_,
  z_] := Module[{T = 1 / sampling, i},
  Table[FullSimplify[hsSecondOrderForm[[i]] /. s -> (2 / T (1 - z^-1) / (1 + z^-1))],
    {i, 1, Length[hsSecondOrderForm]}]
];


getHzImpulseInvariance[hspoles_, coeff_, sampling_, z_] :=
    Module[ {T = 1 / sampling, i, hz, expr, DEBUG = False},
      If[ EvenQ[Length[hspoles]],
        hz = Table[0, {i, Length[hspoles] / 2}];
        Do[
          expr = (z coeff[[i]]) / (z - Exp[hspoles[[i]] T]) + (z  Conjugate@coeff[[i]]) /
              (z - Exp[Conjugate@hspoles[[i]] T]);
          expr = Normal@Chop@FullSimplify[ComplexExpand[expr]];
          expr = ComplexExpand[Numerator[expr]] / ComplexExpand[Denominator[expr]];
          If[ DEBUG,
            Print["expr=", expr ]
          ];
          hz[[i]] = expr;
          If[ DEBUG,
            Print["hz=", hz ]
          ];
          , {i, 1, Length[hz]}
        ],
        hz = Table[0, {i, 1 + ((Length[hspoles] - 1) / 2)}];
        Do[
          expr = (z coeff[[i]]) / (z - Exp[hspoles[[i]] T]) +
              (z  Conjugate@coeff[[i]]) / (z - Exp[Conjugate@hspoles[[i]] T]);
          expr = Normal@Chop@FullSimplify[ComplexExpand[expr]];
          expr = ComplexExpand[Numerator[expr]] / ComplexExpand[Denominator[expr]];
          If[ DEBUG,
            Print["expr=", expr ]
          ];
          hz[[i]] = expr;
          If[ DEBUG,
            Print["hz=", hz ]
          ]
          , {i, 1, Length[hz] - 1}
        ];
        hz[[-1]] = (z coeff[[(Length[hspoles] - 1) / 2 + 1]]) /
            (z - Exp[T hspoles[[(Length[hspoles] - 1) / 2 + 1]]])
      ];
      If[ DEBUG,
        Print["hz=", hz ]
      ];
      Chop@hz
    ];


makehs[poles_, tk_, s_] := Module[{i},
  tk / Product[s - poles[[i]], {i, 1, Length[poles]}]
];


getButterworthQuadraticTerms[hspoles_List, coeff_List, s_] :=
    Module[ {i, hs, expr, DEBUG = False},
      If[ DEBUG,
        Print["enter getButterworthQuadraticTerms"]
      ];
      If[ EvenQ[Length[hspoles]],
        If[ DEBUG,
          Print["EvenQ"]
        ];
        hs = Table[0, {i, Length[hspoles] / 2}];
        Do[
          expr = coeff[[i]] / (s - hspoles[[i]]) + Conjugate@coeff[[i]] /
              (s - Conjugate@hspoles[[i]]);
          expr = Normal@Chop@FullSimplify[ComplexExpand[expr]];
          expr = (Chop@ComplexExpand[Numerator[expr]]) /
              (Chop@ComplexExpand[Denominator[expr]]);
          hs[[i]] = expr
          , {i, 1, Length[hs]}
        ],
        If[ DEBUG,
          Print["OddQ"]
        ];
        hs = Table[0, {i, 1 + ((Length[hspoles] - 1) / 2)}];
        Do[
          If[ DEBUG,
            Print["*********i=", i]
          ];
          expr = coeff[[i]] / (s - hspoles[[i]]) + Conjugate@coeff[[i]] /
              (s - Conjugate@hspoles[[i]]);
          If[ DEBUG,
            Print["expr=", expr]
          ];
          expr = Normal@Chop@FullSimplify[ComplexExpand[expr]];
          If[ DEBUG,
            Print["expr=", expr]
          ];
          expr = Chop@ComplexExpand[Numerator[expr]] /
              (Chop@ComplexExpand[Denominator[expr]]);
          If[ DEBUG,
            Print["expr=", expr]
          ];
          hs[[i]] = expr;
          , {i, 1, Length[hs] - 1}
        ];
        hs[[-1]] = Chop@coeff[[(Length[hspoles] - 1) / 2 + 1]] /
            (s - hspoles[[(Length[hspoles] - 1) / 2 + 1]])
      ];
      If[ DEBUG,
        Print["leaving, hs=", hs]
      ];
      hs
    ];


getButterworthGain[hspoles_List  (* List of stable poles locations *)] := Module[{i},
  Chop[ Product[-hspoles[[i]], {i, 1, Length[hspoles]}]]
];


getButterworthGain[hspoles_List, (* List of stable poles locations *)
  method_?(IntegerQ[#] && Positive[#]&), (* design method, 1=BILINEAR, 2=Impulse Invariance*)
  fs_?(NumberQ[#] && Positive[#]&)        (* sampling frequency in Hz*)
] := Module[{k, i},
  k = Chop[ Product[-hspoles[[i]], {i, 1, Length[hspoles]}]];
  If[method == 1, k, (1 / fs) * k ]
];


getPartialFractionsCoeff[
  poles_List, (*pole locations *)
  gain_                 (*number in the numerator, the H(s) gain *)
] :=
    Module[ {i, s},
      Chop@Table[ Limit[makehs[  Delete[poles, i], gain , s ],
        s -> poles[[i]]], {i, 1, Length[poles]}]
    ];


getHsFromPolesSecondOrder[hspoles_, s_, coeff_] :=
    Module[ {i, hs, expr, DEBUG = False},
      If[ DEBUG,
        Print["enter getHsFromPolesSecondOrder, hspoles=", hspoles, " coeff=", coeff ]
      ];
      If[ EvenQ[Length[hspoles]],
        hs = Table[0, {i, Length[hspoles] / 2}];
        Do[
          expr = coeff[[i]] / (s - hspoles[[i]]) + Conjugate@coeff[[i]] / (s - Conjugate@hspoles[[i]]);
          expr = Chop@Normal@Chop@FullSimplify[ComplexExpand[expr]];
          expr = Chop@ComplexExpand[Numerator[expr]] / ComplexExpand[Denominator[expr]];
          If[ DEBUG,
            Print["expr=", expr ]
          ];
          hs[[i]] = expr;
          If[ DEBUG,
            Print["hs=", hs ]
          ]
          , {i, 1, Length[hs]}],
        hs = Table[0, {i, 1 + ((Length[hspoles] - 1) / 2)}];
        Do[
          expr = coeff[[i]] / (s - hspoles[[i]]) + Conjugate@coeff[[i]] /
              (s - Conjugate@hspoles[[i]]);
          expr = Normal@Chop@FullSimplify[ComplexExpand[expr]];
          expr = ComplexExpand[Numerator[expr]] / ComplexExpand[Denominator[expr]];
          If[ DEBUG,
            Print["expr=", expr ]
          ];
          hs[[i]] = expr;
          If[ DEBUG,
            Print["hs=", hs ]
          ]
          , {i, 1, Length[hs] - 1}
        ];
        hs[[-1]] = coeff[[(Length[hspoles] - 1) / 2 + 1]] /
            (s - hspoles[[(Length[hspoles] - 1) / 2 + 1]]);
      ];
      If[ DEBUG,
        Print["hs=", hs ]
      ];
      hs
    ];


getButterworthStablePoles[
  filterOrder_?(IntegerQ[#] && Positive[#] &), (*filter order, N *)
  omegac_?(NumericQ[#] && Positive[#]&)   (*Radius of butter circle, or cutoff freq.*)
  ] := Module[ {i, poles},
      poles = Table[0, {i, filterOrder}];
      Do[
        poles[[i + 1]] = omegac(Cos[(Pi(1 + 2 i + filterOrder)) / (2 filterOrder)] + I Sin[(Pi(1 + 2 i + filterOrder)) / (2 filterOrder)])
        , 
        {i, 0, Length[poles] - 1}
      ];
      poles
    ];


(*--------------------------------------------*)
(*  decide on phase signa for pulse train     *)
(*--------------------------------------------*)
getPhase[n_, cn_, t0_, period_] :=
    Module[ {phase},
      If[ n == 0,
        Return[0]
      ];
      If[ cn == 0,
        Return[0]
      ];
      If[ t0 == 0,
        Return[If[ cn < 0,
          Sign[n] * Pi,
          0
        ]]
      ];
      phase = -n * (2Pi) / period t0;
      If[ cn >= 0 && Abs[phase] > 0,
        Return[Mod[phase, Sign[phase] * Pi]],
        Return[0]
      ];
      phase = -n * (2Pi) / period t0 + Sign[n] * Pi;
      If[ Abs[phase] > 0,
        Return[Mod[phase, Sign[phase] * Pi]],
        Return[Sign[n] * Pi]
      ];
    ];


convertToLCD[$mat_?(MatrixQ[#, NumericQ] &)] :=
    Module[ {mat = $mat, c},
      mat = Rationalize[mat];
      c = LCM[Sequence @@ Flatten[Denominator[mat]]];
      (Numerator[#] * (c / Denominator[#]) / HoldForm@Evaluate@c) & /@ mat
    ];



(*--------------------------------------------*)
(*  writtenby John Fultz, see                 *)
(* http://mathematica.stackexchange.com/questions/1039/mouseappearance-and-cursor-problems*)
(*--------------------------------------------*)
stripGraphMouseAppearance[x_Graph] :=
    RawBoxes[ToBoxes[x, StandardForm] /.
        TagBox[contents_, MouseAppearanceTag["NetworkGraphics"]] :>
            contents]


(*--------------------------------------------*)
(*  writtenby Hypnotoad                *)
(* http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245*)
(*--------------------------------------------*)
getMaxPadding[p_List] := Map[Max, (BorderDimensions@
    Image[Show[#,
      LabelStyle -> White,
      Background -> White]] & /@ p) ~ Flatten ~ {{3}, {2}}, {2}] + 1;


(*--------------------------------------------*)
(*  writtenby Hypnotoad                *)
(* http://mathematica.stackexchange.com/questions/8241/generating-a-broken-or-snipped-axis-in-listplot/8245#8245*)
(*--------------------------------------------*)
snip[pos_] := Arrowheads[{{Automatic, pos,
  Graphics[{BezierCurve[{{0, -(1 / 2)}, {1 / 2, 0}, {-(1 / 2), 0}, {0, 1 / 2}}]}]}}];


(*--------------------------------------------*)
(*  written by Xerxes from stackexchange   3/3/2013         *)
(*http://mathematica.stackexchange.com/questions/20570/wrapper-for-inexact-numeric-complex-numbers-that-maintains-polar-form*)
(*  *)
(*--------------------------------------------*)
MakeBoxes[polarForm[z_Complex], form_] :=
    With[{r = Abs[z], \[Phi] = Arg[z]},
      RowBox[{If[r == 1, Sequence @@ {}, MakeBoxes[r, form]],
        If[\[Phi] == 0, Sequence @@ {},
          SuperscriptBox[MakeBoxes[E, form],
            RowBox[{MakeBoxes["\[ImaginaryI]", form],
              If[\[Phi] == 1, Sequence @@ {},
                MakeBoxes[\[Phi], form]]}]]]}]];



(*------------------------------*)
checkForSingularity[p_, q_, at_List, x_] := Module[{p0, q0, t, r, r0},
  r = First@Last@Reap@Do[
       If[at[[n]] === Infinity,
        p0 = p /. x -> (1/t);
        p0 = (2 t - p0)/t^2;
        q0 = (q /. x -> (1/t))/t^4;
        r0 = checkForSingularity[p0, q0, {0}, t];
        Sow[Flatten[{Infinity, Rest[Flatten@r0]}]];
        ,
        r0 = {at[[n]], Limit[(x - at[[n]])  p, x -> at[[n]]], 
          Limit[(x - at[[n]])^2 q, x -> at[[n]]]};
        Sow@r0
        ]
       , {n, 1, Length@at}
       ];
  r
  ];
  
End[]
Protect@@Names["nma`*"];
EndPackage[]
