(* :Title: Useful Mathematica functions for Computer Music Theory *)

(* :Author: Craig Stuart Sapp   <craig@ccrma.stanford.edu> *)

(* :Summary:

   This package contains Mathematica functions of use for the courses
   Music 320 and Music 420 at Stanford University in Digital Audio
   Theory and Applications.

*)

(* :Name: SCMTheory.m *)

(* :PackageVersion: 0.00 *)

(* :Version: Mathematica 2.2 *)

(* :Keywords: Computer Music Theory *)

(* :Warnings: none *)

(* :Date: 20 November 1997 *)

(*
 * Changing $Context to SCMTheory` and $ContextPath to {SCMTheory`, System`}:
 *)

Unprotect[Alias];


BeginPackage["SCMTheory`", "Utilities`FilterOptions`"]

(***** Usage Messages ****************************************************)


Alias::usage = 
      "Alias[modSequence, factor] Returns the alias of modSequence by \
       factor.  Length of modSequence divided by factor must be an \
       integer.";

Amplitude::usage = 
      "Amplitude[number] Gives the absolute value of the number which is \
       the length of a complex number.  Slightly smarter than the Abs \
       function.";

CartesianForm::usage = 
      "RectForm[number] Converts a complex number or list of numbers \
       to rectangular form: Re[x] + I Im[x]. Alias: RectForm.";

RectForm::usage = CartesianForm::usage;

CircularPlot::usage = 
      "CircularPlot[modSequence, options] plots the complex number \
       sequence on the domain of a circle.  Magnitude of number is \
       is the distance from the domain circle to the plotted points. \
       Angle of the complex number is the angle of the plotted point \
       on a torroidal ring around the domain circle at the given \
       magnitude radius.  A bandlimited interpolation between \
       points can be added.  Alias is CPlot";
       
CPlot::usage = CircularPlot::usage;

CircularTransformPlot::usage = 
      "CircularTransformPlot[modSequence, options] is similar to the  \
       function CircularPlot, but will plot the DFT of the sequence \
       as points and an interpolation line using the DTFT.";

ComplexSinusoidPlot::usage = 
     "ComplexSinusoidPlot[minTime, maxTime, options] \
      Plots a Complex Sinusoid and projections onto the real and imaginary \
      planes.  Type Options[ComplexSinusoidPlot] to see the default \
      options. Alias is CSPlot";
     
CSPlot::usage = ComplexSinusoidPlot::usage;

Convolution::usage = 
      "Convolution[modSequence1, modSequence2, Spectral->True] \
       Will convolve the two signals.  If Spectral->True, then convolution \
       will be done by multiplying the individual spectral components of the \
       two signals and then taking the inverse DFT. If Spectral->False then \
       convolution will be done in the time domain.  Note that the convolution \
       is circular -- the length of the output sequence is the same as the \
       input sequences.";

Correlation::usage = 
      "Correlation[modSequence1, modSequence2] Multiplies the conjugate of \
       the DFT of the conjugate of sequence1 by the DFT of sequence 2 and \
       returns the inverse DFT of this spectrum.";

Decimate::usage = 
      "Decimate[modSequence, amount] will sample every \
       amount'th value from the modSequence.  modSequence should be divisible \
       by amount+1 or last incomplete multiple of amount in signal will be \
       dropped.  Example: {0,1,2,3} decimated by 1 gives {0,2}.";

DFTAnalyze::usage = 
      "DFTAnalyze[modSequence, options] Treats the input sequence as a \
       signal, and plots the signal and its spectrum.  Uses the function \
       RealImaginaryPlot and MagnitudePhasePlot.";

EvenPart::usage = 
      "EvenPart[function, variable] returns the even part of function with \
       respect to variable.  Even + Odd = original function.";

EvenOddPlot::usage = 
      "EvenOddPlot[function, range, plotOptions...] plots the given function \
       as well as plotting the even and odd parts of the functions over the \
       given range.  Options are the same as for Plot[], plus there is \
       EvenStyle and OddStyle to control the style of the even and odd \
       function lines. Alias: EOPlot";

EOPlot::usage = EvenOddPlot::usage;

EvenOddSinusoidPlot::usage = 
      "EvenOddSinusoidPlot[phaseShift, plotOptions...] plots two cycles \
       of a sine wave with the given phase as a dashed line along with \
       the cosine (even) and sine (odd) functions with no phase that add \
       together to give the phase-shifted sine function.  Same options \
       as EvenOddPlot.  Alias: EOSPlot";

EOSPlot::usage = EvenOddSinusoidPlot::usage;

Flip::usage = 
      "Flip[modSequence] Gives the flip, or reverse, of the modSequence. \
       Example: Flip[{0,1,2,3,4}] = {0,4,3,2,1}";

HermitianQ::usage = 
      "HermitianQ[function] tests whether the function or sequence or \
       2-D matrix is hermitian.";

InnerProduct::usage = 
      "InnerProduct[vector1, vector2] Returns sum of v1 * Conjugate[v2]. \
       Alias: IP.";

IP::usage = InnerProduct::usage;

Interpolate::usage =
      "Interpolate[modSequence, amount, options...] Will stretch the \
       sequence with amount-1 extra samples between the original samples \
       with band-limted values.";

MagnitudePhasePlot::usage = 
      "MagnitudePhasePlot[modSequence, options] plot the magnitude \
       and phase of the input sequenc.  Default display is magnitude \
       in a decibel scale.  Uses the functions MagnitudePlot and \
       PhasePlot.";

MagnitudePlot::usage = 
      "MagnitudePlot[modSequence, options]  plots the Magnitude of a \
       modulo sequence.  Uses the same options as SeqPlot, but adds \
       the option LogScale which plots the y-axis on a decibel scale \
       if set to True (function assumes that sequence contains \
       amplitudes).  Will take continuous functions like SeqPlot, but \
       will only sample the functions and not display them.";

NumericIDFT::usage = 
      "NumericIDFT[modSequence] Returns the numeric Inverse DFT of the \
       sequence.";

NumericDFT::usage =
      "NumericDFT[modSequence] Returns the numeric Discrete Fourier \
       Transform of the sequence.";

OddPart::usage = 
      "OddPart[function, variable] returns the even part of function with \
       respect to variable.  Even + Odd = original function.";

ParabolaFit::usage =
      "ParabolaFit[s1, s2, s3, variable:x] returns a parabolic \
       function of variable that passes through all three specified \
       samples and which is centered about the middle sample \
       ParabolaFit[{x1,y1},{x2,y2},{x3,y3}, variable:x] will return \
       a parabolic function of variable for any three arbitrary points.";

ParabolaPeak::usage = 
      "ParabolaPeak[s1, s2, s3] returns the parabolic extremum \
       point of the parabola uniquely described by the three input \
       samples.  Returns a delta sample value for the location of \
       the peak from the second sample ParabolaPeak[{x1,y1},{x2,y2},{x3,y3}] \
       will give the x-y coordinates of the parabolic maximum (or \
       minimum) of an arbitrary set of three points.";

ParabolaFitPlot::usage = 
      "ParabolaFitPlot[p1,p2,p3, ShowPlot:False, plotOptions___] \
       plots the specified three points or samples and the fitted \
       parabolic function.  If ShowPlot is not False, then the \
       maximum or minimum of the parabola is displayed in red (light \
       gray).  plotOptions are the plotting options similar to Plot's options.";

Phase::usage = 
      "Phase[number] Gives the angle of a complex number.  Slightly smarter \
       than the Arg function.  For computational simplicity, Phase[0] is \
       defined as zero eventhough Arg[0] is undefined.";

PhasePlot::usage = 
      "PhasePlot[modSequence, options] Will plot the phase of a \
       sequence.  Has similar options as SeqPlot, but you cannot \
       use the Repeat option.  See Options[PhasePlot] for default \
       options.";

PolarForm::usage = 
      "Polar[complexNumber] Converts a number or list of numbers to \
       polar form: Abs[x] Exp[I Arg[x]].";

RealImaginaryPlot::usage = 
      "RealImaginaryPlot[modSequence, options...] plots the \
       real part and imaginary part of a sequence.  Uses the \
       two functions RealPlot and ImaginaryPlot.";

ReImPlot::usage = 
      "ReImPlot[function, range, plotOptions...] plots the real part of \
       the given function as a dashed black line and the imaginary \
       part as a solid red line over the given range.  Alias: RIPlot";

RIPlot::usage = ReImPlot::usage;

Repeat::usage = 
      "Repeat[modSequence, repetition] Will repeat modSequence \
       by repetition.  Example: Repeat[{1,2,3}, 2] == {1,2,3,1,2,3}.";

SampleFunction::usage = 
      "SampleFunction[function, {variable, min, max, spacing:1}] Samples \
       a function and returns a list of the samples.";

SeqPlot::usage = 
      "SeqPlot[modSequence, options...] plots a modulo sequence as a discrete \
       set of lines.  SeqPlot[function, range, options...] samples the \
       function, with the possibility of adding the rule Continuous->True \
       to overlay the sampled function.  See Options[SeqPlot] for a list of \
       the default options.";

Shift::usage = 
      "Shift[modSequence, amount] Same as RotateRight[].  Delays a sequence \
       of samples by amount.  Example:  Shift[{0,1,2,3}, 1] == {3,0,1,2}.";

ShiftToCausal::usage = 
      "ShiftToCausal[modSequence] Rotates the sequence by half the \
       length.  If the length is even, then the first element in the \
       sequence is interpreted as the most negative unpaired element.  \
       Reverses exactly the process of ShiftToNonCausal.";

ShiftToNonCausal::usage = 
      "ShiftToNonCausal[modSequence] Rotates the sequence by half the \
       length.  If the length is even, then the unpaired element at \
       length/2 is interpreted as the most negative element.  \
       Reverses exactly the process of ShiftToCausal.";

SignalFunction::usage = 
      "SignalFunction[signalList, variable:n, offset:0] Creates an arbitrary \
       signal where each signal element in the list is an Impulse at time \
       (position+offset), of variable.";

Sinc::usage = 
      "Sinc[x] returns the function Sin[Pi x]/(Pi x).   
       \
       Sinc[length, x] returns the digital sinc function, where \
       length is the separation between an infinite number of \
       shifted sinc functions.";

Stretch::usage = 
      "Stretch[signal, amount] Will stretch a signal by \
       amount.  Example: {1,1,1} stretched by 2 = {1,0,1,0,1,0}.";

SymbolicDFT::usage = 
      "SymbolicDFT[modSequence, variable1:n, variable2:k] Returns a symbolic \
       DFT on the sequence.";

SymbolicIDFT::usage = 
      "SymbolicInverseDFT[modSequence, variable1:k, variable2:n] Returns a \
       symbolic Inverse DFT of the sequence.";

ZeroPad::usage = 
      "ZeroPad[modSequence, numZeros] will zero-pad a sequence with \
       specified number of zeros. Start of signal is assumed to be \
       time 0, and negative time is coming from the other direction, \
       so zeros are added in the middle of the sequence.";


(* Options *)

Options[CartesianForm] := Module[{xxx, yyy}, 
   { Abs[xxx_ + I yyy_]^(-2) -> (xxx^2+yyy^2)^(-1) }]

Options[RectForm] := Options[CartesianForm]

Options[CircularTransformPlot] := {
   EyeSeparation->(5 Degree),
   ZeroPadding->200,
   InterpolationFactor->10,
   Stereo->False,
   ViewPoint->{-2.580, -1.809, 1.143},
   Animation->False,
   Points->True,
   Axis->True,
   Normalize->True,
   ColorPoints->False,
   Bars->True,
   BarColor->True,
   Interpolate->True,
   ColorLine->True,
   Transformation->True,
   ColorAxis->True,
   LineThickness->0.01,
   BarThickness->0.01
};

Options[CTPlot] := Options[CircularTransformPlot]

Options[ComplexSinusoidPlot] := Options[CSPlot]

Options[CSPlot] := {
      (* Physical Variables *)
      Amplitude->1, AngularFrequency -> 2 Pi, 
      Frequency->None, Phase->0, Decay->Infinity,
        
      (* Plotting Options: *)
      ShowSinusoid->True, ShowComplex->True, 
      ShowReal->True, ShowImaginary->True,
      ShowConnectors->True, ShowProjection->True, 
      PlotPoints->Automatic,
   
      RealConnector->True, ImaginaryConnector->True, 
      ComplexConnector->True, ComplexPosition->-1.1,
      
      Boxed->True, SinusoidBoxed->True, 
      ProjectionBoxed->True,
      SinusoidBoxStyle->Thickness[0.0],
      ProjectionBoxStyle->RGBColor[0.7,0.7,0.7],
      ProjectionJoinStyle->Dashing[{0.01, 0.01}],
      BoxScale->Automatic, ComplexBoxed -> True,  
      RealBoxed->True,
      ImaginaryBoxed->True, SinusoidBoxed->Automatic,
   
      RealStyle->RGBColor[0, 0.5, 0], 
      ImaginaryStyle->RGBColor[0, 0, 0.5],
      ComplexStyle->RGBColor[0.5, 0, 0], 
      SinusoidStyle->Thickness[0.005],
      ConnectStyle->{RGBColor[0.3, 0.3, 0.3], Dashing[{0.01, 0.01}]},
   
      RealOffset->{0, 1}, ImaginaryOffset->{1, 0}, 
      ComplexOffset->{-1,-1},
    
      System`ViewPoint->{-10, -24, -6}
   }


Options[Convolution] := {Spectral->True}


Options[DFTAnalyze] := {SignalOptions->{Causal->True}, 
         SpectrumOptions->{Causal->False}, 
         SignalPlot->RealImaginaryPlot, SpectrumPlot->MagnitudePhasePlot}

Options[EOPlotOdd] := {OddStyle->xRed}
Options[EOPlotEven] := {EvenStyle->xGreen}
Options[EvenOddPlot] := Flatten[{Options[EOPlotOdd], Options[EOPlotEven]}]
Options[EOPlot] := Options[EvenOddPlot]

Options[EvenOddSinusoidPlot] := Options[EvenOddPlot]
Options[EOSPlot] := Options[EvenOddPlot]

Options[ImaginaryPlot] := {Frame->True, Axes ->{Automatic, None}, 
         FrameTicks->{Automatic, Automatic, None, None},
         Interpolate->True, InterpStyle->Thickness[0.0], 
         FrameLabel->{"", "", 
            FontForm["Imaginary Part", {"Times-Roman", 18}], ""}, 
         Stretch->1, Lines -> False}

Options[MagnitudePhasePlot] := {LeftOptions->{}, RightOptions->{}}

Options[MagnitudePlot] := {LogScale->True, Frame->True, Axes->False,
         FrameTicks->{Automatic, Automatic, None, None},
         Interpolate->True, InterpStyle->Thickness[0.0], 
         FrameLabel->{"","",FontForm["Magnitude", {"Times-Roman", 18}], ""},
         Stretch->1, Lines->False}

Options[ParabolaFitPlot] := {Peak->True, PeakStyle->{RGBColor[1, 0, 0],
         PointSize[0.02]}, PointStyle->PointSize[0.02]}

Options[PhasePlot] := {Frame->True, Axes->{True, None}, PlotRange->All, 
         Interpolate->True, InterpStyle->Thickness[0.0],
         Stretch->1, Wrap->True, Frame->True, Lines->False,
         PointStyle->PointSize[0.025], BreakStyle->Dashing[{0.005,0.03}],
         FrameLabel->{"","",FontForm["Phase", {"Times-Roman", 18}], ""} }

Options[RealImaginaryPlot] := {Axes->{True, None}, LeftOptions->{}, 
         RightOptions->{}}

Options[RealPlot] := {Frame->True, Axes ->{Automatic, None}, 
         FrameTicks->{Automatic, Automatic, None, None},
         Interpolate->True, Stretch->1, Lines -> False, 
         InterpStyle->Thickness[0.0], PlotRange->All,
         FrameLabel->{"", "", FontForm["Real Part", {"Times-Roman", 18}], ""} }

Options[SeqPlot] = {Points->True, Lines->True, 
         PointStyle->PointSize[0.03], LineStyle->Thickness[0.012],   
         Repeat->1, Causal->True, TickPeriod->Automatic, Continuous->False,
         PlotRange->All, Stretch->1, PadFactor->12, 
         InterpStyle->Thickness[0.02], Interpolate->False, InterpGraphic->Line}
  



(***** Start of Function Definitions *************************************)


Begin["`Private`"]


(***** Alias **************************************************************)

Alias[list_?System`ListQ, factor_] := Module[{len = Length[list]},
   Map[Apply[Plus, #]&,Transpose[Partition[list, len/factor]]]
 /; IntegerQ[len/factor] ]



(***** Amplitude **********************************************************)

Amplitude[x_] := Module[ {rform = RectForm[x]},
   Simplify[Sqrt[Im[rform]^2 + Re[rform]^2]] ]

AppendTo[Attributes[Simplify], Listable];
(* Print["Listable added to Simplify[]'s attributes."] *)
AppendTo[Attributes[Amplitude], Listable];



(***** CartesianForm *****************************************************)

CartesianForm[x_, options___] := Module[ 
   { both = Together[ComplexExpand[ReleaseHold[x]]] },
   ComplexExpand[Simplify[Re[both]] + 
      Factor[Simplify[I Im[both]]] ] /. {options} /. Options[CartesianForm]
]

RectForm[parameters___] := CartesianForm[parameters]

AppendTo[Attributes[CartesianForm], Listable];
AppendTo[Attributes[RectForm], Listable];



(***** CircularTransformPlot *********************************************)

CircularTransformPlot[signalList_?System`ListQ, opts___] := Module[
   {
      r = 3.0,
      x1, y1, z1,
      x2, y2, z2,

      signal = N[signalList],
      length = Length[signalList],
      normalizationFactor,
      discreteStartingPoints, discreteEndingPoints, discreteValues,
      plot1, plot2, fixPlot,
      discretePointsG,
      discreteBarsG,
      interpolateG,
      unitCircleG,
      options
   },
   options = Flatten[{opts, Options[CircularTransformPlot]}];

   zeroPadding = ZeroPadding /. options;
   eyeSeparation = EyeSeparation /. options;
   viewPoint = ViewPoint /. options;
   interpolationFactor = InterpolationFactor /. options;
   view = View /. options;

   If[Transformation /. options,
      discreteValues = NumericDFT[N[signal]],
      discreteValues = N[signal]];

   (* Normalize the discreteValues: *)
   If[Normalize /. options,
      normalizationFactor = Max[Abs[discreteValues]];
      discreteValues /= normalizationFactor
   ];

   (*** Create the graphics for discrete items ***)
   discreteStartingPoints = Table[ N[{Cos[2Pi i/length], Sin[2Pi i/length],
      0}], {i, 0, length-1}];
   discreteEndingPoints = TorroidalSurfaceMap[discreteValues];
   If[Points /. options,
      discretePointsG = {PointSize[0.02], Map[Point,
         discreteEndingPoints]},
      discretePointsG = {}
   ];

   If[Bars /. options,
      Block[ {tempLines = Map[Line, Transpose[{discreteStartingPoints,
         discreteEndingPoints}]]},
      discreteBarsG = {Thickness[BarThickness /. options],
         If[BarColor /. options,
            Transpose[{Map[Hue[Chop[Phase[#1]/(2. Pi)]]&, discreteValues],
               tempLines}], {RGBColor[0,0,0], tempLines} ]}
      ]
      ,
      discreteBarsG = {}
   ];



   (*** Make the interpolated spectral line: ***)
   If[Interpolate /. options,
      Block[ {tempLine, interpolatedHues, fragmentedLine, mappedPoints,
         interpolatedLength},
         If[Transformation /. options,
            interpolatedLength = length + zeroPadding,
            interpolatedLength = length * interpolationFactor
         ];

         If[Transformation /. options,
            interpolatedValues = NumericDFT[ZeroPad[signal, zeroPadding]],
            interpolatedValues = Interpolate[signal, interpolationFactor]
         ];
         If[Normalize /. options, interpolatedValues /= normalizationFactor];
  
         mappedPoints = TorroidalSurfaceMap[interpolatedValues];
         tempLine = Join[mappedPoints, {First[mappedPoints]}];
         If[ColorLine /. options,
            interpolatedHues = Map[Hue[Chop[Phase[#1]/(2. Pi)]]&,
               interpolatedValues];
            fragmentedLine = Map[Line, Partition[Drop[Drop
                  [Flatten[Map[{#,#}&, tempLine], 1], -1], 1], 2]];
            interpolateG = {Thickness[LineThickness /. options],
               Transpose[{interpolatedHues, fragmentedLine}]}
            , (* else *)
            interpolateG = {Thickness[LineThickness /. options],
               RGBColor[0.01,0.04,0.1], Line[tempLine]}
         ]
      ],
      interpolateG = {} (* not interpolating *)
   ];


   (*** Graphics for the axis ***)
   If[Axis /. options,
      Block[ {$DisplayFunction = Identity},
         If[ColorAxis /. options,
            unitCircleG = ParametricPlot3D[{Cos[t], Sin[t], 0, Hue[t/(2Pi)]},
               {t,0,2Pi}, Axes->False],
            unitCircleG = ParametricPlot3D[{Cos[t], Sin[t], 0},
               {t,0,2Pi}, Axes->False]
         ]
      ],
      unitCircleG = {} (* no axis *)
   ];


   (* Check for Fixed Plots for animations *)
   If[Animation /. options,
      fixPlot = {PointSize[0.0], RGBColor[1,1,1], Point[{-2,-2,-1}],
         Point[{2,2,1}]},
      fixPlot = {}
   ];


   (* Determine the ViewPoint: *)
   If[ Length[viewPoint] === 2,
         x1 = r Cos[viewPoint[[1]]];
         y1 = r Sin[viewPoint[[1]]];
         z1 = viewPoint[[2]];
         x2 = r Cos[viewPoint[[1]] + eyeSeparation];
         y2 = r Sin[viewPoint[[1]] + eyeSeparation];
         z2 = z1;
      ,
         x1 = viewPoint[[1]];
         y1 = viewPoint[[2]];
         z1 = viewPoint[[3]];
         x2 = viewPoint[[1]] + eyeSeparation;
         y2 = viewPoint[[2]] + eyeSeparation;
         z2 = z1;
   ];

   (* Create the 3D plot *)
   plot1 = Graphics3D[{discretePointsG, discreteBarsG, interpolateG, fixPlot}];

   (* finally, display the plots *)
   If[Stereo /. options,
      Block[{$DisplayFunction = Identity},
         plot2 = {Show[Graphics3D[{discretePointsG, discreteBarsG,
            interpolateG, fixPlot}], unitCircleG, ViewPoint->{x1, y1, z1},
            PlotRange->All, FilterOptions[Plot3D, options]],
         Show[Graphics3D[{discretePointsG, discreteBarsG,
            interpolateG, fixPlot}], unitCircleG, ViewPoint->{x2, y2, z2},
            PlotRange->All, FilterOptions[Plot3D, options]]}
      ];
      Show[GraphicsArray[plot2]]
   ,
      (* Mono view *)
      Show[plot1, unitCircleG, ViewPoint->{x1, y1, z1}, PlotRange->All,
         FilterOptions[Plot3D, options]]
   ] 
]

TorroidalSurfaceMap[list_?ListQ] := Module[ {length = Length[list]},
   MapIndexed[ N[{
         Cos[2 Pi (First[#2]-1)/length] (1 + Abs[#1] Cos[Pi/2+Phase[#1]]),
         Sin[2 Pi (First[#2]-1)/length] (1 + Abs[#1] Cos[Pi/2+Phase[#1]]),
         Abs[#1] Sin[Pi/2+Phase[#1]]
      }]&,
      list]
]

CTPlot[parameters__] := CircularTransformPlot[parameters]

CircularPlot[signal_, options___] :=
   CircularTransformPlot[signal, options,
      Transformation->False]

CPlot[parameters__] :=  CircularPlot[parameters]



(***** ComplexSinusoidPlot ***********************************************)

CSin[t_, A_, omega_, phi_, tau_] := A E^(-t/tau) E^(I(t omega + phi))

ComplexSinusoidPlot[min_, max_, options___] := Module[
   {
      amp, omega, phi, tau,
      xvalues, yvalues, zvalues,
      maxReal, maxImaginary,
      sinusoidPoints, 
      realPoints, 
      imaginaryPoints, 
      complexPoints, 
      scaleRe, scaleIm,
      sinusoidPlot = {}, realPlot = {}, imaginaryPlot = {}, complexPlot = {},
      sinusoidBox = {}, reBox = {}, imBox = {}, compBox = {},
      reConnect = {}, imConnect = {}, compConnect = {},
      sinusoidStyle = {}, realStyle = {}, imaginaryStyle = {},complexStyle = {},
      range = max - min, fstep,
      viewPoint
   }, 

   (* Setup defaults *)
   amp = Amplitude /. {options} /. Options[CSPlot];
   omega = AngularFrequency /. {options} /. Options[CSPlot];
   frequency = Frequency /. {options} /. Options[CSPlot];
   If[frequency =!= None, omega = 2 Pi frequency];
   phi = Phase /. {options} /. Options[CSPlot];
   tau = Decay /. {options} /. Options[CSPlot];
   sinusoidStyle = SinusoidStyle /. {options} /. Options[CSPlot];
   realStyle = RealStyle /. {options} /. Options[CSPlot];
   imaginaryStyle = ImaginaryStyle /. {options} /. Options[CSPlot];
   complexStyle = ComplexStyle /. {options} /. Options[CSPlot];
   connectStyle = ConnectStyle /. {options} /. Options[CSPlot];
   viewPoint = ViewPoint /. {options} /. Options[CSPlot];
   viewPoint = Flatten[viewPoint];
   boxScale = BoxScale /. {options} /. Options[CSPlot];
   sinusoidBoxStyle = SinusoidBoxStyle /. {options} /. Options[CSPlot];
   boxed = Boxed /. {options} /. Options[CSPlot];
   sinusoidBoxed = SinusoidBoxed /. {options} /. Options[CSPlot];
   showSinusoid = ShowSinusoid /. {options} /. Options[CSPlot];
   showReal = ShowReal /. {options} /. Options[CSPlot];
   showImaginary = ShowImaginary /. {options} /. Options[CSPlot];
   showComplex = ShowComplex /. {options} /. Options[CSPlot];
   showConnectors = ShowConnectors /. {options} /. Options[CSPlot];
   showProjection = ShowProjection /. {options} /. Options[CSPlot];
   complexPosition = ComplexPosition /. {options} /. Options[CSPlot];
   plotPoints = PlotPoints /. {options} /. Options[CSPlot];
   projectionBoxStyle = ProjectionBoxStyle /. {options} /. Options[CSPlot];
   projectionBoxed = ProjectionBoxed /. {options} /. Options[CSPlot];


   (* Generate the plotting points *)
   Which[
      NumberQ[plotPoints], fstep = range / plotPoints,
      True, fstep = range / 200.0
   ];
   xvalues = SampleFunction[Re[CSin[t, amp, omega, phi, tau]], 
      {t, min, max, fstep}];
   maxReal = Max[Abs[xvalues]];
   yvalues = SampleFunction[Im[CSin[t, amp, omega, phi, tau]], 
      {t, min, max, fstep}];
   maxImag = Max[Abs[yvalues]];
   zvalues = Range[min, max, fstep];
   zeros = zvalues * 0;

   (* calculate the scaling factor for the boxes *)
   Which[
      (* If just one number, scale the size of box *)
      NumberQ[boxScale], 
         scaleRe = maxReal * boxScale; scaleIm = maxImag * boxScale,
      (* If two numbers use as absolute values of box lengths (animations) *)
      ListQ[boxScale] && Length[boxScale] === 2, 
         scaleRe = boxScale[[1]]; scaleIm = boxScale[[2]],
      (* If not scaling fit sinusioid exactly in box *)
      True, scaleRe = maxReal; scaleIm = maxImag
   ];

   (* Draw the Complex Sinusoid *)
   If[showSinusoid,
      sinusoidPoints = Transpose[{zvalues, -xvalues, yvalues}];
      sinusoidPlot = Flatten[{ SinusoidStyle /. Options[CSPlot],
            sinusoidStyle, Line[sinusoidPoints] }];
   ];

   (* Draw the imaginary part of the sinusoid *)
   If[showImaginary && showProjection,
      imaginaryPoints = Transpose[{zvalues,zeros-scaleRe, yvalues+2.5*scaleIm}];
      imaginaryPlot = Flatten[{ ImaginaryStyle /. Options[CSPlot], 
         imaginaryStyle, Line[imaginaryPoints] }];
   ];

   (* Draw the real part of the sinusoid *)
   If[showReal && showProjection,
      realPoints = Transpose[{zvalues,-xvalues+4*scaleRe, zeros-scaleIm}];
    (*  realPoints = Transpose[{zvalues,-xvalues+4*scaleRe, zeros-scaleRe}]; *)
      realPlot = Flatten[{ RealStyle /. Options[CSPlot], 
         realStyle, Line[realPoints] }];
   ];

   (* Draw the complex part of the sinusoid *)
   complexPosition = complexPosition * (scaleRe + scaleIm) / 2.0;
   If[showComplex && showProjection,
      complexPoints = Transpose[{min+zeros+complexPosition, 
         -xvalues, yvalues}];
      complexPlot = Flatten[{ ComplexStyle /. Options[CSPlot], 
         complexStyle, Line[complexPoints] }];
   ];

  
   (* Draw the Complex Sinusoid box *)
   If[boxed && sinusoidBoxed, 
      sinusoidBox = Flatten[{sinusoidBoxStyle,
         Line[{{min, -scaleRe, -scaleIm}, 
            {min, scaleRe, -scaleIm}, {min, scaleRe, scaleIm}, 
            {min, -scaleRe, scaleIm}, {min, -scaleRe, -scaleIm}, 
            {max, -scaleRe, -scaleIm}}],
         Line[{{max, scaleRe, scaleIm}, 
            {max, -scaleRe, scaleIm}, {max, -scaleRe, -scaleIm}, 
            {max, scaleRe, -scaleIm}, {max, scaleRe, scaleIm}, 
            {min, scaleRe, scaleIm}}],
         Line[{{max, -scaleRe, scaleIm}, {min, -scaleRe, scaleIm}}], 
         Line[{{max, scaleRe, -scaleIm}, {min, scaleRe, -scaleIm}}]
       }];
   ];

  
   If[projectionBoxed && boxed && showProjection, 
      (* Draw the Imaginary box *)
      imBox =  Flatten[{ProjectionBoxStyle /. Options[CSPlot], 
         projectionBoxStyle,
         Line[{ {min, -1*scaleRe, 1.5*scaleIm},
                {max, -1*scaleRe, 1.5*scaleIm},
                {max, -1*scaleRe, 3.5*scaleIm},
                {min, -1*scaleRe, 3.5*scaleIm},
                {min, -1*scaleRe, 1.5*scaleIm} }]  }];
      (* Draw the Real box *)
      reBox =  Flatten[{ProjectionBoxStyle /. Options[CSPlot], 
         projectionBoxStyle,
         Line[{ {min, 3*scaleRe, -1*scaleIm},
                {max, 3*scaleRe, -1*scaleIm},
                {max, 5*scaleRe, -1*scaleIm},
                {min, 5*scaleRe, -1*scaleIm},
                {min, 3*scaleRe, -1*scaleIm} }]  }];
      (* Draw the Complex box *)
      compBox =  Flatten[{ProjectionBoxStyle /. Options[CSPlot], 
         projectionBoxStyle,
         Line[{ {complexPosition+min, scaleRe, scaleIm},
                {complexPosition+min, -scaleRe, scaleIm},
                {complexPosition+min, -scaleRe, -scaleIm},
                {complexPosition+min, scaleRe,-scaleIm},
                {complexPosition+min, scaleRe, scaleIm} }]  }];
   ];


   If[showConnectors && boxed && projectionBoxed && showProjection, 
      reConnect = Flatten[{ConnectStyle /. Options[CSPlot], connectStyle,
         Line[{{min, scaleRe, -scaleIm}, {min, 3*scaleRe, -scaleIm}}],
         Line[{{max, scaleRe, -scaleIm}, {max, 3*scaleRe, -scaleIm}}] }];
      imConnect = Flatten[{ConnectStyle /. Options[CSPlot], connectStyle,
         Line[{{min, -scaleRe, scaleIm}, {min, -scaleRe, 1.5*scaleIm}}],
         Line[{{max, -scaleRe, scaleIm}, {max, -scaleRe, 1.5*scaleIm}}] }];
      compConnect = Flatten[{ConnectStyle /. Options[CSPlot], connectStyle,
         Line[{{min, scaleRe, scaleIm}, 
            {complexPosition+min, scaleRe, scaleIm}}],
         Line[{{min, -scaleRe, scaleIm}, 
            {complexPosition+min, -scaleRe, scaleIm}}],
         Line[{{min, -scaleRe, -scaleIm}, 
            {complexPosition+min, -scaleRe, -scaleIm}}],
         Line[{{min, scaleRe, -scaleIm}, 
            {complexPosition+min, scaleRe, -scaleIm}}] }];
   ];


   Show[Graphics3D[Flatten[{sinusoidPlot, sinusoidBox, realPlot, reBox,
      imaginaryPlot, imBox, complexPlot, compBox, reConnect,
      imConnect, compConnect }]], 
      Boxed->False, PlotRange->All, ViewPoint->viewPoint] 
]


CSPlot[parameters___] := ComplexSinusoidPlot[parameters]




(***** Convolution *******************************************************)

Convolution[modSequence1_, modSequence2_, options___] :=
   If[Spectral /. {options} /. Options[Convolution],
      NumericIDFT[NumericDFT[modSequence1] NumericDFT[modSequence2]]
   ,  Apply[Plus, Transpose[Table[RotateRight[modSequence1, i] modSequence2,
         {i, 0, Length[modSequence2]-1}]]]
]



(***** Correlation *******************************************************)

Correlation[signalList1_?System`ListQ, signalList2_?System`ListQ] :=
   NumericIDFT[Conjugate[NumericDFT[Conjugate[signalList1]]] *
      NumericDFT[signalList2]]



(***** Decimate **********************************************************)

Decimate[modSequence_List, amount_?Positive]:=
   First[Transpose[Partition[modSequence, amount]]]



(***** DFTAnalyze ********************************************************)

DFTAnalyze[modSequence_?VectorQ, options___] := Module[
   {
      plot1, plot2, 
      transform = NumericDFT[modSequence],
      signalOptions = SignalOptions /. {options} /. Options[DFTAnalyze], 
      spectrumOptions = SpectrumOptions /. {options} /. Options[DFTAnalyze] 
   },
   Block[{$DisplayFunction = Identity},
      plot1 = (SignalPlot /. {options} /. Options[DFTAnalyze])[modSequence,
                 Evaluate[Apply[Sequence, signalOptions]],
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence, Options[DFTAnalyze]]]  
              ];
      plot1 = Show[plot1, Frame->True,
                 FrameLabel->{"",FontForm["signal", {"Times-Roman", 18}], "",""}
              ];
      plot2 = (SpectrumPlot /. {options} /. Options[DFTAnalyze])[transform,
                 Evaluate[Apply[Sequence, spectrumOptions]],
                 Evaluate[Sequence[options]],      
                 Evaluate[Apply[Sequence, Options[DFTAnalyze]]] 
              ];  
      plot2 = Show[plot2, Frame->True,
                 FrameLabel->{"",FontForm["spectrum",{"Times-Roman", 18}],"",""}
              ];
   ];
   Show[plot1];
   Show[plot2]
]



(***** EvenPart **********************************************************)

EvenPart[modSequence_List] := 1/2 (modSequence + Flip[modSequence])

EvenPart[function_, variable_] := 
   1/2 (function + (function /.  variable-> -variable))



(***** EvenOddPlot *******************************************************)

EvenOddPlot[function_, {variable_, min_, max_},  options___] := Module[
   { 
     plainPlot, evenPlot, oddPlot, 
     evenStyle = PlotStyle->(EvenStyle /. {options} /. Options[EOPlotEven]),
     oddStyle = PlotStyle->(OddStyle /. {options} /. Options[EOPlotOdd])
   },
   Block[{$DisplayFunction=Identity},
      plainPlot = Plot[Evaluate[function], {variable, min, max},
         Evaluate[FilterOptions[Plot, options]],
         Evaluate[PlotStyle->xDashing]
         ];
      evenPlot = Plot[Evaluate[EvenPart[function, variable]], 
         {variable, min, max},
         Evaluate[evenStyle]
         ];
      oddPlot = Plot[Evaluate[OddPart[function, variable]], 
         {variable, min, max},
         Evaluate[oddStyle]
         ]; 
   ];
   Show[plainPlot, evenPlot, oddPlot]
]

EOPlot[parameters___] := EvenOddPlot[parameters]



(***** EvenOddSinusoidPlot ***********************************************)


EvenOddSinusoidPlot[phase_, options___] := Module[ {x},
   EvenOddPlot[Sin[2Pi x + phase], {x, -1, 1}, 
      options,
      Frame->True,
      DefaultFont->{"Symbol", 18},
      FrameTicks->{{{-1, "-2\[Pi]"}, {-0.5, "-\[Pi]"}, 
         {0, "0"}, {0.5, "\[Pi]"}, {1, "2\[Pi]"}}, {-1,1}, None, None}]
]

EOSPlot[parameters___] := EvenOddSinusoidPlot[parameters] 



(***** Flip **************************************************************)

Flip[modSequence_?ListQ] := RotateRight[Reverse[modSequence]]



(***** HermitianQ ********************************************************)


HermitianQ[list_?VectorQ] := Module[
   {
      re = Apply[And,Map[(# == 0)&, Re[OddPart[list]]]],
      im = Apply[And,Map[(# == 0)&, Im[EvenPart[list]]]]
   },
   re && im
]


HermitianQ[list_?MatrixQ] := Module[
   {i,j,dim=Dimensions[list]},
   If[Length[dim] == 2 && dim[[1]] == dim[[2]], 
   For[i=1, i<=Length[list], i++,
      For[j=i, j>=1, j--,
         If[list[[i]][[j]] != Conjugate[list[[j]][[i]]], Return[False]]
      ]
   ];
   Return[True],
   False
   ]
]


HermitianQ[function_, variable_] := Module[{},
   ((function /. variable->-variable) /. {
      a_ b_ -> Conjugate[a] b}) === function]



(***** ImaginaryPlot: see RealPlot ***************************************)

(***** InnerProduct ******************************************************)

InnerProduct[v1_?VectorQ, v2_?VectorQ] := Apply[Plus, v1 Conjugate[v2]]

IP[parameters___] := InnerProduct[parameters]



(***** Interpolate *******************************************************)

Interpolate[modSequence_List, amount_?Positive] := Module[
   { analysis = NumericDFT[modSequence], length = Length[modSequence] },
   amount NumericIDFT[ZeroPad[analysis, length * (amount-1)]] 
]


(***** MagnitudePlot *****************************************************)

MagnitudePhasePlot[modSequence_?VectorQ, options___] := Module[
   {
      plot1, plot2,
      leftOptions = LeftOptions /. {options} /. Options[MagnitudePhasePlot],
      rightOptions = RightOptions /. {options} /. Options[MagnitudePhasePlot]
   },
   Block[{$DisplayFunction = Identity},
      plot1 = MagnitudePlot[modSequence,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence, leftOptions]],
                 Evaluate[Apply[Sequence, Options[MagnitudePhasePlot]]] 
               ];
      plot2 = PhasePlot[modSequence,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence, rightOptions]],
                 Evaluate[Apply[Sequence, Options[MagnitudePhasePlot]]] 
               ];
   ];
   Show[GraphicsArray[{plot1, plot2}], GraphicsSpacing->0]
]



(***** MagnitudePlot *****************************************************)

MagnitudePlot[modSeq_?VectorQ, options___] := Module[
   {
      interp, interpSeq, plot1, plot2, newSeq,
      length = Length[modSeq],
      logScale = LogScale /. {options} /. Options[MagnitudePlot],
      interpolate = Interpolate /. {options} /. Options[MagnitudePlot],
      stretch = Stretch /. {options} /. Options[MagnitudePlot]
   },
   Which[ 
      length > 256, interp = 1, (* no need to interp. for line if many points *)
      length > 128, interp = 2,
      length > 64,  interp = 4,
      length > 32,  interp = 8,
      True,         interp = 16
   ];
   interpSeq = Interpolate[modSeq, interp];
   
   (* plot1 is for discrete samples, plot2 is for interpolation if requested *)
   
   Block[{$DisplayFunction = Identity}, 
   If[logScale,
   (* Log scale plots *)
   newSeq = Chop[Map[20 Log[10, #]&, N[Abs[modSeq]]]];
   newSeq = Map[If[# == -Infinity, -1000, #]&, newSeq]; (* replace -Infinity *)
   plot1 = SeqPlot[newSeq,
              Interpolate->False,
              Evaluate[Sequence[options]],
              Evaluate[Apply[Sequence,Options[MagnitudePlot]]],
              Evaluate[If[logScale,
                 PlotRange->Automatic, PlotRange->All]]
           ];
   newInterpSeq = Chop[Map[20 Log[10, #]&, N[Abs[interpSeq]]]];
   newInterpSeq = Map[If[# == -Infinity, -1000, #]&, newInterpSeq]; 
   plot2 = If[interpolate,
              SeqPlot[newInterpSeq,
                 Interpolate->True, PadFactor->1,
                 Lines->False, Points->False, Stretch->stretch/interp,
                 Ticks->None,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence,Options[MagnitudePlot]]]
              ]
           ,
              {}
           ];
   , (* Linear scale plot *)
   plot1 = SeqPlot[Chop[N[Abs[modSeq]]],
              Interpolate->False,
              Evaluate[Sequence[options]],
              Evaluate[Apply[Sequence,Options[MagnitudePlot]]],
              Evaluate[If[logScale,
                 PlotRange->Automatic, PlotRange->All]]
           ];
           
   plot2 = If[interpolate,
              SeqPlot[Chop[N[Abs[interpSeq]]],
                 Interpolate->True, PadFactor->1, Ticks->None,
                 Lines->False, Points->False, Stretch->stretch/interp,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence,Options[MagnitudePlot]]]
              ]
           ,
              {}
           ];
    ]; (* end if logscale *)
    ]; (* end of no display block *)
 
   Show[{plot1, plot2}]

]
      
MagnitudePlot[function_, {var_, min_, max_, step_:1}, options___] :=
   MagnitudePlot[SampleFunction[function, {var, min, max, step}], options]



(***** NumericIDFT *******************************************************)

NumericIDFT[modSequence_?ListQ] := 
   If[AllZeroQ[modSequence], modSequence,
      Chop[N[InverseFourier[Flip[modSequence]] / Sqrt[Length[modSequence]]]]]



(***** NumericDFT ********************************************************)

NumericDFT[modSequence_?ListQ] := 
   If[AllZeroQ[modSequence], modSequence, 
      Chop[N[Sqrt[Length[modSequence]] Fourier[Flip[modSequence]]]]]

AllZeroQ[list_?VectorQ] := N[list[[1]]] == 0 && Apply[Equal, N[list]]



(***** OddPart ***********************************************************)

OddPart[modSequence_?ListQ] := 1/2 (modSequence - Flip[modSequence])


OddPart[function_, variable_]  := 
   1/2 (function - (function /.  variable-> -variable));



(***** ParabolaFit ********************************************************)

ParabolaFit[{x1_:-1,y1_}, {x2_:0,y2_}, {x3_:1,y3_}, var_:Global`x]:=
   Module[ {a,b,p,solutions},
      solutions = Solve[{a(x1-p)^2+b==y1,a(x2-p)^2+b==y2,a(x3-p)^2+b==y3},
         {a,b,p}];
      a (var - p)^2 + b /. solutions
   ]

ParabolaFit[y1_?NumberQ, y2_?NumberQ, y3_?NumberQ] :=
   ParabolaFit[{y1}, {y2}, {y3}]



(***** ParabolaPeak *******************************************************)

ParabolaPeak[s1_?NumberQ, s2_?NumberQ, s3_?NumberQ] :=
   1/2 (s1 - s3) / (s1 - 2 s2 + s3)

ParabolaPeak[{x1_,y1_}, {x2_,y2_}, {x3_,y3_}] :=
   Module[ {
      rules = Solve[{a(x1-p)^2+b==y1, a(x2-p)^2+b==y2, a(x3-p)^2+b==y3},
         {a,b,p}]
      },
      Flatten[{p,b} /. rules]
   ]



(***** Phase *************************************************************)

Phase[x_] := Module[ {rform = RectForm[x]},
   If[rform === 0 || rform === 0., 0,
      ArcTan[Re[rform], Im[rform]]] //. Options[PolarForm] 
]

AppendTo[Attributes[Phase], Listable]



(***** PhasePlot ***********************************************************)

PhasePlot[modSeq_?VectorQ, options___] := Module[
   {
      interp, interpSeq,
      plot1, plot2, branch,
      length = Length[modSeq],
      stretch = Stretch /. {options} /. Options[PhasePlot],
      causal  = Causal  /. {options} /. Options[PhasePlot] /. Options[SeqPlot],
      interpolate = Interpolate /. {options} /. Options[PhasePlot],
      breakStyle = Flatten[{BreakStyle} /. {options} /. Options[PhasePlot]] 
   },
   Which[ 
      length > 256, interp = 2, 
      length > 128, interp = 4,
      length > 64,  interp = 8,
      length > 32,  interp = 16,
      True,         interp = 64
   ];
   
   (* plot1 is for discrete samples, plot2 is for interpolation if requested *)
   
   Block[{$DisplayFunction = Identity},  
   If[Wrap /. {options} /. Options[PhasePlot],
    
     (* Wrapped phase plot *)
      plot1 = SeqPlot[Chop[Phase[N[modSeq]]],
              Interpolate->False, Repeat->1,
              Evaluate[Sequence[options]],
              Evaluate[Apply[Sequence,Options[PhasePlot]]],
              FrameTicks->{Automatic, { {-Pi, tpi["-\[Pi]"]}, 
                 0, {Pi, tpi["\[Pi]"]} }, None, None }
              ]; 
      If[interpolate, 
         interpSeq = Chop[Phase[N[Interpolate[modSeq, interp]]]];
         list = If[causal,
             Transpose[{Range[0, length*interp-1] stretch/interp, interpSeq}],
             Transpose[{Range[-Floor[length*interp-1]/2,
                Floor[length*interp-1]/2] * stretch/interp, 
                Shift[interpSeq, Floor[interp length/2]]}]
         ];
         plot2 = Graphics[Map[If[(N[Abs[(#[[1,2]])-(#[[2,2]])]] > N[0.9 Pi]) ||
            N[#[[2,2]]] == 0 || N[#[[1,2]]] == 0, 
            {Evaluate[Apply[Sequence,breakStyle]],Line[#]}, Line[#]]&, 
            Partition[Rest[Rest[RotateRight[Flatten[Map[{#, #}&, list], 
            1]]]],2]]]
      ,
         plot2 = {}
      ];
      
   , (* unwrapped phase plot *)  
   
   interpSeq = Chop[N[Phase[Interpolate[modSeq, interp]]]];
   
   list = If[causal,
             Transpose[{Range[0, length*interp-1] stretch/interp, interpSeq}],
             Transpose[{Range[-Floor[length*interp-1]/2, 
                Floor[length*interp-1]/2] * stretch/interp, 
                Shift[interpSeq, Floor[interp length/2]]}]
          ];

   (* unwrapping the phase *)
   newModSeq = If[causal,
                  N[Phase[modSeq]],
                  N[Phase[Shift[modSeq, Floor[length/2]]]]
               ];
               
   For[i=2; branch=0; last=list[[1,2]], i < length interp + 1, i++,
      Which[
         list[[i,2]]-last > N[0.75 2Pi], branch--,
         last-list[[i,2]] > N[0.75 2Pi], branch++
      ];
      If[ IntegerQ[(i-1)/interp+1], newModSeq[[(i-1)/interp+1]] += branch 2Pi ];
      last = list[[i,2]];
      list[[i,2]] = list[[i, 2]] + branch 2Pi;
   ];
   
   newModSeq = If[causal,
                  newModSeq,
                  Shift[newModSeq, -Floor[length/2]]
               ];
   plot1 = SeqPlot[newModSeq,
           Interpolate->False, Repeat->1,
           Evaluate[Sequence[options]],
           Axes->None,  
           Evaluate[Apply[Sequence,Options[PhasePlot]]],
           FrameTicks->{Automatic, {{-Pi, tpi["-\[Pi]"]}, 
           0, {Pi, tpi["\[Pi]"]}, {2Pi, tpi["2\[Pi]"]}, {-2Pi, tpi["-2\[Pi]"]},
           {3Pi, tpi["3\[Pi]"]}, {-3Pi, tpi["-3\[Pi]"]}, {4Pi, tpi["4\[Pi]"]},
           {-4Pi, tpi["-4\[Pi]"]}, {5Pi, tpi["5\[Pi]"]}, {-5Pi, tpi["-5\[Pi]"]},
           {6Pi, tpi["6\[Pi]"]}, {-6Pi, tpi["-6\[Pi]"]}, {7Pi, tpi["7\[Pi]"]},
           {-7Pi, tpi["-7\[Pi]"]}, {8Pi, tpi["8\[Pi]"]}, {-8Pi, tpi["-8\[Pi]"]},
           {9Pi, tpi["9\[Pi]"]}, {-9Pi, tpi["-9\[Pi]"]}, {10Pi, tpi["10\[Pi]"]},
           {-10Pi, tpi["-10\[Pi]"]}, {11Pi, tpi["11\[Pi]"]}  
           }, None, None}
        ]; 
   If[interpolate, 
      plot2 = Map[If[(N[Abs[(#[[1,2]])-(#[[2,2]])]] > N[0.9 Pi]) || 
            N[Mod[#[[2,2]],2Pi]] == 0 || N[Mod[#[[1,2]], 2Pi]] == 0,
            {Evaluate[Apply[Sequence,breakStyle]], Line[#]}, Line[#]]&, 
            Partition[Rest[Rest[RotateRight[Flatten[Map[{#, #}&, list], 
            1]]]],2]];
      plot2 = Graphics[plot2];
   ,
      plot2 = {}
   ];

   ]; (* end if wrapping *)
       
   (* making sure to plot the full range from -Pi to Pi in branch 0 *)
   plot3 = Show[Graphics[{PointSize[0], RGBColor[1,1,1], 
              Point[{0,Pi*1.05}], Point[{0,-Pi*1.05}]}]];

   ];  (* end block *)
   
   Show[{plot1, plot2, plot1, plot3}]
]
      
      
PhasePlot[function_, {var_, min_, max_, step_:1}, options___] :=
   PhasePlot[SampleFunction[function, {var, min, max, step}], options]


tpi[string_] := FontForm[string, {"Symbol", 18}] 



(***** ParabolaFitPlot *****************************************************)

ParabolaFitPlot[{x1_:-1,y1_}, {x2_:0,y2_}, {x3_:1,y3_}, options___] := Module[
   {fit, plot, mixplot, peak, x, peakStyle, showPeak, pointStyle},
   peakStyle = PeakStyle /. {options} /. Options[ParabolaFitPlot];
   showPeak = Peak/. {options} /. Options[ParabolaFitPlot];
   pointStyle = PointStyle /. {options} /. Options[ParabolaFitPlot];
   peak = {}; 
   If[showPeak,
      peak = Flatten[{PeakStyle /. Options[ParabolaFitPlot], peakStyle,
         Point[ ParabolaPeak[{x1,y1},{x2,y2},{x3,y3}] ]}]  
   ];
   fit = ParabolaFit[{x1,y1},{x2,y2},{x3,y3}, x];
   Block[{$DisplayFunction = Identity},
      plot = 
        Plot[fit, {x, x1-1, x3+1}, Evaluate[FilterOptions[Plot, options]]];
      mixplot = Show[{plot,
         Graphics[{peak, PointStyle /. Options[ParabolaFitPlot], pointStyle,
            Map[Point,{{x1,y1},{x2,y2},{x3,y3}}]
         }]
      }];
   ];
   Show[mixplot]
]
			

ParabolaFitPlot[y1_?NumberQ, y2_?NumberQ, y3_?NumberQ, options___] := 
   ParabolaFitPlot[{y1}, {y2}, {y3}, options]



(***** PolarForm *********************************************************)

PolarForm[x_, options___] := Module[
   {
      exponent =
         If[(Re[x] =!= 0) && (Re[x] =!= 0.0),
            Exp[If[Head[Arg[x]]===Real,
               HoldForm[Evaluate[ I N[Arg[x]/Pi] Pi]],
               I Arg[x] ]],
            Exp[HoldForm[Evaluate[I Arg[x]]]] ],
      rules, length
   },
   rules = {a->ExpandAll[Abs[x]],b->exponent}//.{options}//.Options[PolarForm];
   length = Simplify[a /. rules];
   If[length === 1 || length === 1.,
      HoldForm[b] /. rules,
      If[exponent === 1 || exponent === 1.,
         HoldForm[a] /. rules,
         If[exponent === -1 || exponent === -1.,
            HoldForm[-a] /. rules,
            HoldForm[a b] /. rules]]]
]

(* next line makes the function work with lists of numbers: *)
AppendTo[Attributes[PolarForm], Listable]

(* Here are some extra rules to help simplify certain expressions: *)

Options[PolarForm] := Module[{x, y}, 
   {
      ArcTan[-x_] -> -ArcTan[x], (* This line must be first for Phase[] *)
      Arg[x_] -> ArcTan[Im[x] / Re[x]],
      Abs[x_] -> Sqrt[Re[x]^2 + Im[x]^2],
      (x_ Cos[y_]^2 + x_ Sin[y_]^2) -> x,
      ArcTan[Tan[x_]] -> x,
      (* ArcTan[-x_] -> -ArcTan[x], *)
      ArcTan[2-Sqrt[3]] -> Pi/12,
      ArcTan[1/Sqrt[3]] -> Pi/6,
      ArcTan[1] -> Pi/4,
      ArcTan[Sqrt[3]] -> Pi/3,
      ArcTan[2+Sqrt[3]] -> 5 Pi/12,
      ArcTan[Infinity] -> Pi/2,
      ArcTan[-2-Sqrt[3]] -> 7 Pi/12,
      ArcTan[-Sqrt[3]] -> 2 Pi/3,
      ArcTan[Sqrt[2]] -> Pi/8,
      ArcTan[1+Sqrt[2]] -> 3 Pi/8,
      ArcTan[1-Sqrt[2]] -> -Pi/8,
      Exp[I ArcTan[Cot[x_]]] -> -Exp[-I(x+Pi/2)],
      Im[x_+I y_]->y, Re[x_+I y_]->x,
      Cos[y_]^2 +  Sin[y_]^2 -> 1,
      x_ Cos[y_]^2 + x_ Sin[y_]^2 -> x,
      Sqrt[(x_)^2] -> x
   }];



(***** RealImaginaryPlot *************************************************)

RealImaginaryPlot[modSequence_?VectorQ, options___] := Module[
   {
      plot1, plot2,
      leftOptions = LeftOptions /. {options} /. Options[RealImaginaryPlot],
      rightOptions = RightOptions /. {options} /. Options[RealImaginaryPlot]
   },
   Block[{$DisplayFunction = Identity},
      plot1 = RealPlot[modSequence,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence, leftOptions]],
                 Evaluate[Apply[Sequence, Options[RealImaginaryPlot]]] 
               ];
      plot2 = ImaginaryPlot[modSequence,
                 Evaluate[Sequence[options]],
                 Evaluate[Apply[Sequence, rightOptions]],
                 Evaluate[Apply[Sequence, Options[RealImaginaryPlot]]] 
               ];
   ];
   Show[GraphicsArray[{plot1, plot2}], GraphicsSpacing->0]
]



(***** RealPlot *** ImaginaryPlot ****************************************)

RealPlot[modSeq_?VectorQ, options___] := Module[
   {
      interp, interpSeq, newSeq, plot1, plot2,
      length = Length[modSeq],
      interpolate = Interpolate /. {options} /. Options[RealPlot],
      stretch = Stretch /. {options} /. Options[RealPlot]
   },
   Which[
      length > 256, interp = 1, (* no need to interp. for line if many points *)
      length > 128, interp = 2,
      length > 64,  interp = 4,
      length > 32,  interp = 8,
      True,         interp = 16
   ];
   interpSeq = Interpolate[modSeq, interp];
  
   (* plot1 is for discrete samples, plot2 is for interpolation if requested *)
  
   Block[{$DisplayFunction = Identity},
      plot1 = SeqPlot[Chop[N[Re[modSeq]]],
              Interpolate->False,
              Evaluate[Sequence[options]],
              Evaluate[Apply[Sequence,Options[RealPlot]]]
           ];
      plot2 = If[interpolate,
                 SeqPlot[Chop[N[Re[interpSeq]]],
                    Interpolate->True, PadFactor->1, Ticks->None,
                    Lines->False, Points->False, Stretch->stretch/interp,
                    Evaluate[Sequence[options]],
                    Evaluate[Apply[Sequence,Options[RealPlot]]]
                 ]
              ,
                 {}
              ];
    ]; (* end of no display block *)

   Show[{plot1, plot2}]

]
     
RealPlot[function_, {var_, min_, max_, step_:1}, options___] :=
   RealPlot[SampleFunction[function, {var, min, max, step}], options]

ImaginaryPlot[modSequence_?VectorQ, options___] :=
   RealPlot[ I modSequence, options, 
      Evaluate[Apply[Sequence, Options[ImaginaryPlot]]]]
   
ImaginaryPlot[function_, {var_, min_, max_, step_:1}, options___] :=
   ImaginaryPlot[SampleFunction[function, {var, min, max, step}], options]



(***** ReImPlot **********************************************************)

ReImPlot[function_, {variable_, min_, max_}, options___] := Module[
   {rePlot, imPlot},
   Block[{$DisplayFunction=Identity},
      rePlot = Plot[Evaluate[Re[function]], {variable, min, max}, 
         options,
         PlotStyle->Dashing[{0.03,0.01}]
         ];
      imPlot = Plot[Evaluate[Im[function]], {variable, min, max},
         PlotStyle->xRed
         ];
   ];
   Show[rePlot, imPlot]
]

RIPlot[parameters___] := ReImPlot[parameters]



(***** Repeat ************************************************************)

Repeat[modSequence_List, amount_?IntegerQ] := Module[ { output={} },
   For[i=0, i<amount, i++, output = Join[output, modSequence]];
   output
]



(***** SampleFunction ****************************************************)

SampleFunction[function_, {variable_, min_, max_, step_:1}] := 
   Table[function, Evaluate[{variable, min, max, step}]]



(***** SeqPlot ***********************************************************)

SeqPlot[list_List, options___] := Module[
   {
     dummyPlot, length, lineStyle, lines, maxRange, minRange, newList, 
     pointLocations, pointStyle, points, range, repeat, stretch, xticks, 
     interpLine
   },
   minRange = Min[list] * 1.05;       (* Adding blank space above and below *)
   maxRange = Max[list] * 1.05;
   repeat = Repeat /. {options} /. Options[SeqPlot];
   stretch = Stretch /. {options} /. Options[SeqPlot];
   newList = Repeat[list, repeat];
   length = Length[newList];
   If[Causal /. {options} /. Options[SeqPlot], 
      range = Range[0,length-1];
   ,  newList = RotateRight[newList, Ceiling[length/2]-1];
      range = Range[Floor[1-length/2], Floor[length/2]];
      If[ Mod[length, 2] == 0, 
         newList = RotateRight[newList];
         range = RotateRight[range];
         range[[1]] = -range[[1]];
      ];
   ];
   range = stretch range;
   pointLocations = Transpose[{range, newList}];
   pointStyle = Flatten[{{PointStyle /. Options[SeqPlot]}, 
      {PointStyle /. {options}/.  Options[SeqPlot]}}];
   points = If[Points /. {options} /. Options[SeqPlot],
               Join[Flatten[{pointStyle}], Map[Point, pointLocations]]
            ,  {} ];
   lineStyle = LineStyle /. {options} /. Options[SeqPlot];
   lines = If[Lines /. {options} /. Options[SeqPlot],
              Join[Flatten[{lineStyle}], Map[Line[{#, {#[[1]], 0}}]& , 
                 pointLocations]]
           ,  {} ];
   xticks = MakeTicks[TickPeriod/.{options}/.Options[SeqPlot], length, stretch];
   interpLine = If[Interpolate /. {options} /. Options[SeqPlot],
                   MakeInterpolationLine[list, options], {} ];
   Block[{$DisplayFunction = Identity},
      dummyPlot = Plot[0,{x,0,0.0000001}, 
         Evaluate[FilterOptions[Plot, options]], 
         Evaluate[If[Frame /. {options} /. Options[Plot],
               FrameTicks->{xticks, Automatic, None, None}
            ,  Ticks->{xticks, Automatic}]],
         DefaultFont->{"Times-Roman", 16}, 
         PlotRange->All,
         PlotStyle->{RGBColor[1,1,1]}  
         ];
   ];
   
   Show[dummyPlot, Graphics[{interpLine, {PointSize[0], RGBColor[1,1,1], 
      Point[{0.0000001,maxRange}],
      Point[{0.0000001,minRange}]}, 
      lines, points}]]
]


SeqPlot[function_, {var_, min_, max_, step_:1}, options___] := Module[
   {sequence, plot1, plot2, len, stretch},
   stretch = Stretch /. {options} /. Options[SeqPlot];
   sequence = SampleFunction[function, {var, min, max, step}];
   Block[{$DisplayFunction = Identity},
      plot1 = SeqPlot[sequence, options];
   ];
   If[Continuous /. {options} /. Options[SeqPlot],
      len = Length[sequence];
      Block[{$DisplayFunction = Identity},
         plot2 = SeqContinuous[function, {var, min, max, step}, 
                    len, stretch, options];
      ];
      Show[plot1, plot2, plot1]   (* ordering 1-2-1 places cont. line in back *)
   ,  Show[plot1]
   ]
]


MakeInterpolationLine[list_, options___] := Module[
   {
      padfactor = PadFactor /. {options} /. Options[SeqPlot],
      repeat = Repeat /. {options} /. Options[SeqPlot],
      causal = Causal /. {options} /. Options[SeqPlot],
      stretch = Stretch /. {options} /. Options[SeqPlot],
      lineStyle = Flatten[{{InterpStyle /. Options[SeqPlot]},
         {InterpStyle /. {options}/.  Options[SeqPlot]}}],
      interpGraphic = InterpGraphic /. {options} /. Options[SeqPlot],

      length, interp, points, range, line

   },
   length = Length[list] padfactor repeat;
   interpolation = Repeat[Re[Interpolate[list, padfactor]], repeat];
   If[causal, 
      range = Range[0,length-1];
   ,  interpolation = RotateRight[interpolation, Ceiling[length/2]];
      range = Range[0, length-1] - Ceiling[length/2]+1;
      range = RotateRight[range];
      range[[1]] = -range[[1]];
   ];
   range = N[stretch range / padfactor];
   points = Transpose[{range, interpolation}];

   If[(interpGraphic === Line),
      line = Join[lineStyle, {Line[points]}],
      line = Join[lineStyle, Map[Point, points]]
   ];
   line
]


SeqContinuous[f_, {var_, min_, max_, step_:1}, len_, stretch_, options___] := Module[ 
   { 
      repeat = Repeat /. {options} /. Options[SeqPlot],
      function = f /. var->(var step/stretch + min), 
      newlen = stretch len
   },
   If[Causal /. {options} /. Options[SeqPlot],
      Table[
         Plot[Evaluate[function /. var->(var-shift)],
            Evaluate[{var, shift, shift + newlen}],
            Evaluate[FilterOptions[Plot, options]]],
         {shift, 0, (repeat-1) newlen, newlen}
      ]      
      ,
      If[ Mod[len repeat, 2] != 0,
         { 
            Table[
               Plot[Evaluate[function /. var->(var-shift)],
                  Evaluate[{var, shift, shift+newlen/2}],
                  Evaluate[FilterOptions[Plot, options]]],
               {shift, -Floor[(repeat-1)/2] newlen, (repeat-1)newlen/2, newlen}
            ]      
            ,  
            Table[
               Plot[Evaluate[function /. var->(var-shift)],
                  Evaluate[{var, shift+newlen/2, shift+newlen/2+newlen/2}],
                  Evaluate[FilterOptions[Plot, options]]], 
               {shift, -(repeat-1)/2 newlen - newlen, (repeat-1)newlen/2-newlen, newlen}
            ]  
         }
         ,
         {
            Table[
               Plot[Evaluate[function /. var->(var-shift)],
                  Evaluate[{var, shift, shift+ newlen/2}],
                  Evaluate[FilterOptions[Plot, options]]],
               {shift, -Floor[(repeat)/2] newlen, (repeat-1)newlen/2, newlen}
            ]      
         ,
            If[ Mod[repeat, 2] == 0,
               Table[
                  Plot[Evaluate[function /. var->(var-shift)],
                     Evaluate[{var, shift+newlen, shift+newlen/2}],
                     Evaluate[FilterOptions[Plot, options]]],
                  {shift, -Floor[(repeat)/2] newlen, (repeat-1)newlen/2, newlen}
               ]
            ,
               Table[
                  Plot[Evaluate[function /. var->(var-shift)],
                     Evaluate[{var, shift+newlen, shift+newlen/2}],
                     Evaluate[FilterOptions[Plot, options]]],
                  {shift, -Ceiling[(repeat)/2] newlen, (repeat-2)newlen/2, newlen}
               ]
            ]
         }
      ]
   ]
]


MakeTicks[period_?NumberQ, length_, stretch_] := Module[ {},
   Flatten[{Table[tick, {tick, 0, -Ceiling[length/2] stretch, -period stretch}],
            Table[tick, {tick, 0, (length-1) stretch, period stretch}]}]
]

MakeTicks[automatic_Symbol, length_, stretch_] := Module[
   { nperiod = If[stretch > 1, 1, Round[1/stretch]] },
   Which[
      length <= 13, MakeTicks[1, length, stretch],
      length <= 25, MakeTicks[2, length, stretch],
      length <= 51, MakeTicks[5, length, stretch],
      length <= 100, MakeTicks[10, length, stretch],
      length <= 180, MakeTicks[20, length, stretch],
      True, MakeTicks[Round[length/10/10]*10, length, stretch]
   ]
]




(***** Shift *************************************************************)

Shift[modSequence_List, amount_?IntegerQ] := RotateRight[modSequence, amount]



(***** ShiftToCausal *****************************************************)

ShiftToCausal[modSeq_?VectorQ , length_:-Infinity] := Module[ {len},
   If[length == -Infinity, len = Length[modSeq], len = length];
   RotateLeft[modSeq, Floor[len/2]]
]



(***** ShiftToNonCausal **************************************************)

ShiftToNonCausal[modSeq_?VectorQ , length_:-Infinity] := Module[ {len},
   If[length == -Infinity, len = Length[modSeq], len = length];
   RotateRight[modSeq, Floor[len/2]]
]



(***** SignalFunction ****************************************************)

SignalFunction[list_?System`ListQ, variable_:Global`n, offset_:0] :=
   GSequenceToFunction[list, variable, offset]



(***** Stretch ***********************************************************)

Sinc[x_] := If[x == 0, 1.0, N[Sin[Pi x]/(Pi x)]]


Sinc[len_, x_] := Module[ {},
   If[OddQ[len], 
      N[ Sin[Pi x]/(len Sin[Pi x / len])]
   ,
      N[Sin[Pi x] Cos[Pi x / len] / (len Sin[Pi x /len])]
   ]
]


DigitalSinc[len_, x_] := Sinc[len, x]



(***** Stretch ***********************************************************)

Stretch[modSequence_List, amount_?IntegerQ] := Module[ 
   { length = Length[modSequence] },
   If[amount <= 1, modSequence,
      Flatten[Transpose[Join[{modSequence},
         Partition[Table[0,{i,length*(amount-1)}], length]]]]]]



(***** SymbolicDFT *******************************************************)

SymbolicDFT[signalList_, variable1_:Global`n, variable2_:Global`k] := Module[
   {
      length = Length[signalList],
      output = {},
      signal, transform
   },
   signal = SignalFunction[signalList, variable1];
   transform = GDFTransform[signal, length, variable1, variable2];
   rules = {SignalProcessing`Support`SigProc`L->length, variable2:>i};

   For[i=0, i<length, i++,
      AppendTo[output, First[transform /. rules]]];

   output
]


(***** SymbolicIDFT ******************************************************)

SymbolicIDFT[spectralList_, variable1_:Global`k, variable2_:Global`n] :=
   SymbolicDFT[Flip[spectralList/Length[spectralList]], variable1, variable2]



(***** ZeroPad ***********************************************************)

ZeroPad[modSequence_List, numZeros_?IntegerQ] := Pad[modSequence, numZeros, 0]

Pad[modSequence_List, num_?IntegerQ, value_:0] := Module[ 
   { length = Length[modSequence], firstPart, secondPart },
   If[num < 1, 
      modSequence
   ,  firstPart = Take[modSequence, Ceiling[length/2]];
      secondPart = Take[modSequence, -Floor[length/2]];
      Join[firstPart, Table[value,{i,num}], secondPart]
   ]
]





(***** Start of Private Definitions **************************************)

xRed = RGBColor[0.8, 0.3, 0.3];
xGreen = RGBColor[0.3, 0.8, 0.3];
xBlue = RGBColor[0.3, 0.3, 0.8];
xDashing = Dashing[{0.02, 0.01}];

FilterOptions[command_Symbol, opts___] := Module[
   {keywords = First /@ Options[command]},
   Sequence @@ Select[ {opts}, MemberQ[keywords, First[#]]& ]
]


GSequenceToFunction[parameters___] :=
   SignalProcessing`Support`SigProc`SequenceToFunction[parameters]

GDFTransform[parameters___] :=
   SignalProcessing`Digital`DFT`DFTransform[parameters]

GInvDFTransform[parameters___] :=
   SignalProcessing`Digital`DFT`InvDFTransform[parameters]


(***** End of Function Definitions ***************************************)

End[] (* SCMTheory`Private` *)


(* Proctect all of the defined functions (when package is finalized):

Protect[
   EOPlot,
   EOSPlot,
   EvenOddPlot,
   EvenOddSinusoidPlot,
   EvenPart,
   OddPart,
   ZeroPad
]

*)

EndPackage[]

(* :Examples:
 
   CartesianForm[Sqrt[2] E^(I Pi/4)]
   Convolution[{1,1,1,1,0,0,0,0}, {1,0,0,0,0,1,1,1}]
   Correlation,
   Decimate[{0,1,2,3,4,5,6,7,8,9}, 2]
   EOPlot[E^(-(x+1)^2), {x, -4, 4}, Ticks->None];
   EOSPlot[Pi/3];
   EvenPart[Cos[x + 1/2]]
   EvenPart[{0,1,2,3}]
   Flip[{0,1,2,3,4}]
   HermitianQ[{1,I,I,-I,-I}]
   HermitianQ[E^(I x), x]
   HermitianQ[{{1,I,2I},{-I,2,-3I},{-2I,3I,3}}]
   Interpolate[{1,0,0,0,0}, 2]
   NumericDFT[{0,1,2,3}]
   NumericIDFT[{7,-2-3I,-3,-2+3I}]
   OddPart[Cos[x + 1/2]]
   OddPart[{0,1,2,3}]
   Phase[I+1]
   PolarForm[I+1]
   RectForm[Sqrt[2] E^(I Pi/4)]
   Repeat[{1,2,3},4]
   RIPlot[E^(I x), {x, -Pi, Pi}, Ticks->None];
   SampleFunction[Sin[2Pi x/8]]
   SeqPlot[{0,1,2,3,2,1}];
   Shift[{0,1,2,3}, 2]
   SignalFunction[{0,1,2,3}]
   Stretch[{1,1,1}, 3]
   SymbolicDFT[{1,1,1,0,0,0}]
   SymbolicIDFT[{7,-2-3I,-3,-2+3I}]
   ZeroPad[{3,2,1,2}, 3]

 *)

