(* :Title: Useful Mathematica functions for generating windows *)

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

(* :Summary:

   This package contains Mathematica functions of use for generating
   analysis windows for Music 420 at CCRMA, Stanford University.

*)

(* :Name: SoundFile.m *)

(* :PackageVersion: 0.00 *)

(* :Version: Mathematica 3.0 *)

(* :Keywords: Sound Soundfile *)

(* :Warnings: none *)

(* :Date: 22 April 1997 *)

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

Needs["Utilities`FilterOptions`"]

(* BeginPackage["SCMT`Windows`"] *)

(***** Usage Statements **************************************************)

Blackman::usage = 
      "BlackmanW[length, options] returns the Blackman window.";

BlackmanW::usage = Blackman::usage;
BlackmanF::usage = Blackman::usage;

CosineW::usage = 
      "CosineW[length, options] return a cosine^Power window.  Default
       option: Power->2.";

RectangularF::usage = 
      "RectangularF[var] the rectangle function.";

RectangularW::usage = 
      "RectangularF[legnth, options] returns a rectangular window.";

RectangleW::usage = RectangularW::usage;

TriangleF::usage = 
      "TriangleF[var] the triangle function.";

TriangleW::usage = 
      "TriangleF[legnth, options] returns a triangle window.";

Window::usage = 
      "Window[length, options] returns a window.";


(***** Options ***********************************************************)

Options[BlackmanHarris] := {A0->0.35875, A1->0.48829, A2->0.14128, A3->0.01168} 
Options[Blackman] := {A0->7938/18608, A1->9240/18608, A2->1430/18608} 
Options[Cauchy] := {Alpha->2.0}
Options[Cosine] := {Power->1}
Options[DolphChebyshev] := {Alpha->3.0}
Options[Gaussian] := {Alpha->2.0}
Options[Hamming] := {Alpha->25/46}
Options[HannPoisson] := Options[Poisson]
Options[Kaiser] := {Alpha->2.0}
Options[Poisson] := {Alpha->2.0}
Options[Turkey] := {Taper->0.5}
Options[Window] := {Type->Rectangular, Causal->True}

Options[WindowTransformPlot] := {Length->51, Normalize->True, Title->None}
Options[WindowPlot] := {Title->None}


$Windows = {Blackman, BlackmanHarris, Bohman, Cauchy, Cosine,
            Gaussian, Hamming, Hann, HannPoisson, Kaiser, 
            Poisson, Rectangle, Reisz, Riemann, Triangle, 
            Turkey, VallePoussin 
           };



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

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


(***** Blackman ***********************************************************)

BlackmanF[x_, options___] := Module[
   {
      a0 = A0 /. {options} /. Options[Blackman],
      a1 = A1 /. {options} /. Options[Blackman],
      a2 = A2 /. {options} /. Options[Blackman]
   },
   N[If[Abs[x]<=1, a0 + a1 Cos[Pi x] + a2 Cos[2Pi x], 0]]
]

BlackmanW[length_, options___] := Window[length, Type->Blackman, options]



(***** BlackmanHarris *****************************************************)

BlackmanHarrisF[x_, options___] := Module[
   {
      a0 = A0 /. {options} /. Options[BlackmanHarris],
      a1 = A1 /. {options} /.  Options[BlackmanHarris],
      a2 = A2 /. {options} /.  Options[BlackmanHarris],
      a3 = A3 /. {options}/. Options[BlackmanHarris]
   },
   N[If[Abs[x]<=1, a0-a1 Cos[Pi (x-1)]+a2 Cos[2Pi(x-1)]-a3 Cos[3Pi(x-1)], 0]]
]

BlackmanHarrisW[length_, options___] := Window[length, Type->BlackmanHarris, 
   options]



(***** Bohman *************************************************************)

BohmanF[x_, options___] := N[If[Abs[x]<=1, 
   (1.0-Abs[x])Cos[Pi x]+1/Pi Sin[Pi Abs[x]], 0]
]

BohmanW[length_, options___] := Window[length, Type->Bohman, options]



(***** Cauchy *************************************************************)

CauchyF[x_, options___] := Module[
   {
      alpha = Alpha /. {options} /. Options[Cauchy]
   },
   N[If[Abs[x]<=1, 1/(1+(alpha x)^2), 0]]
]

CauchyW[length_, options___] := Window[length, Type->Cauchy, options]



(***** Cosine *************************************************************)

CosineF[x_, options___] := Module[
   {
      power = Power /. {options} /. Options[Cosine]
   },
   N[If[Abs[x]<=1, Cos[Pi x/2]^power, 0]]
]
   
CosineW[length_, options___] := Window[length, Type->Cosine, options]

CosineLobeW[length_, options___] := CosineW[length, Power->1, options]
RaisedCosineW[length, options___] := CosineW[length, Power->2, options]



(***** DolphChebyshev *****************************************************)

DolphChebyshevW[length_, options___] := Module[
   {
      alpha = Alpha /. {options} /.
      Options[DolphChebyshev],
      beta,
      spectrum,
      len = length-1 
   },
   beta[len_, alpha_] := Cosh[1/len ArcCosh[10^alpha]];
   spectrum = Table[Chop[N[
      (-1)^k Cos[len  ArcCos[B[len, alpha] Cos[Pi k/len]]]/
      Cosh[len ArcCosh[B[len, alpha]]]]], 
      {k, 0, len-1}];
   signal = NumericIDFT[spectrum];
   signal /= Max[signal];
   output = AppendTo[signal, First[signal]];
   If[Causal /. {options} /. Options[Window],
      output,
      RotateRight[ ShiftToNonCausal[output], 1]
   ]
]



(***** Gaussian ***********************************************************)

GaussianF[x_, options___] := Module[
   {
      alpha = Alpha /. {options} /.  Options[Gaussian]
   },
   N[If[Abs[x]<=1, Exp[-(alpha x)^2/2], 0]]
]

GaussianW[length_, options___] := Window[length, Type->Gaussian, options]



(***** Hamming ************************************************************)

HammingF[x_, options___] := Module[
   {
      alpha = Alpha /. {options} /. Options[Hamming]
   },
   N[If[Abs[x]<=1, alpha+(1-alpha) Cos[x Pi], 0]]
]

HammingW[length_, options___] := Window[length, Type->Hamming, options]



(***** Hann ***************************************************************)

HannF[x_, options___] := N[If[Abs[x]<=1, 1/2 + 1/2 Cos[x Pi], 0]]

HannW[length_, options___] := Window[length, Type->Hann, options]



(***** HannPoisson ********************************************************)

HannPoissonF[x_, options___] := HannF[x, options] PoissonF[x, options]

HannPoissonW[length_, options___] := Window[length, Type->HannPoisson, options]



(***** Kaiser *************************************************************)

KaiserF[x_, options___] := Module[
   {
      alpha = Alpha /. {options} /. Options[Kaiser]
   },
   N[If[Abs[x]<=1, BesselI[0,Pi alpha Sqrt[1.0-x^2]]/BesselI[0,Pi alpha], 0]]
]

KaiserW[length_, options___] := Window[length, Type->Kaiser, options]



(***** Poisson ************************************************************)

PoissonF[x_, options___] := Module[
   {
      alpha = Alpha /. {options} /.  Options[Poisson]
   },
   N[If[Abs[x]<=1, Exp[-alpha Abs[x]], 0]]
]

PoissonW[length_, options___] := Window[length, Type->Poisson, options]



(***** Rectangular ********************************************************)

RectangleF[x_, options___] := N[If[Abs[x]<=1, 1, 0]]

RectangleW[length_, options___] := Window[length, Type->Rectangle, options]

RectangularW[parameters___] := RectangleW[parameters]
DirichletW[parameters___] := RectangleW[parameters]



(***** Reisz **************************************************************)

ReiszF[x_, options___] := N[If[Abs[x]<=1, 1-Abs[x]^2, 0]]

ReiszW[length_, options___] := Window[length, Type->Reisz, options]



(***** Riemann ************************************************************)

RiemannF[x_, options___] := N[If[Abs[x]<=1, Sinc[x], 0]]

RiemannW[length_, options___] := Window[length, Type->Riemann, options]



(***** Triangle ***********************************************************)

TriangleF[x_, options___] := N[If[Abs[x]<=1, 1-Abs[x], 0]]

TriangleW[length_, options___] := Window[length, Type->Triangle, options]

BartletW[parameters___] := TriangleW[parameters]
FejerW[parameters___] := TriangleW[parameters]



(***** Turkey *************************************************************)

TurkeyF[x_, options___] := Module[
   {
      alpha = 1.0 - Taper /. {options} /.  Options[Turkey]
   },
   N[Which[
      Abs[x] <= alpha,  1, 
      Abs[x] <= 1,      1/2 (1+Cos[Pi(Abs[x]-alpha)/(1-alpha)]),
      True,             0
   ]]
]

TurkeyW[length_, options___] := Window[length, Type->Turkey, options]



(***** VallePoussin *******************************************************)

VallePoussinF[x_, options___] := N[Which[
   Abs[x] <= 1/2,  1-6x^2(1-Abs[x]),
   Abs[x] <= 1,    2(1-Abs[x])^3,
   True,           0
]]

VallePoussinW[length_, options___] := Window[length, Type->VallePoussin, 
   options]



(***** Window *************************************************************)

Window[length_, options___] := Module[
   {
      name = Symbol[ToString[Type /. {options} /. Options[Window]] <> "F"],
      output
   },
   output = N[SampleFunction[name[x, options], {x, -1, 1 , 2.0/(length-1)}]];
   If[Causal /. {options}/. Options[Window],
      output,
      ShiftToCausal[output, length]
   ]
]

(***** WindowPlot *********************************************************)

WindowPlot[windowType_Symbol, options___] := Module[
   {
      titleTest = Title /. {options} /.
      Options[WindowPlot],
      name = Symbol[ToString[windowType] <> "F"],
      plotOptions,
      titlex
   },
   If[titleTest === None,
      titlex = ToString[windowType] <> " Window",
      titlex = ToString[titleTest];
   ];
   plotOptions[title_] :=Sequence[ {x, -1.2, 1.2},
      Frame->True, 
      AspectRatio->Automatic,
      PlotRange->{{-1.2,1.2}, {-0.1,1.1}}, 
      Axes->False,
      FrameTicks->{Automatic,{0,0.25, 0.5, 0.75, 1.0}, None, None},
      DefaultFont->{"Times-Roman", 14},
      FrameLabel->{"","", FontForm[title, {"Times-Roman",
      18}], ""}
   ];
   Plot[name[x, options], Evaluate[plotOptions[Evaluate[titlex]]], 
      Evaluate[FilterOptions[Plot,options]]]
]



(***** WindowTransformPlot ************************************************)

WindowTransformPlot[windowType_Symbol, options___] := Module[
   {
      length = Length /. {options} /.  Options[WindowTransformPlot],
      titleTest = Title /. {options} /.  Options[WindowTransformPlot],
      name = Symbol[ToString[windowType] <> "W"],
      transform
   },
   If[titleTest == None,
      title = ToString[windowType]<> " Window Transform",
      title = titleTest
   ];
   transform = NumericDFT[name[length, Causal->False, options]]/length;
   If[Normalize /. {options} /.
      Options[WindowTransformPlot],
      transform /= Abs[Max[Re[transform]]]];
   MagnitudePlot[  
      transform, 
      Evaluate[Sequence[options]],
      Lines->False,
      FrameLabel->{"","dB", FontForm[title,
      {"Times-Roman", 18}], ""},
      PointStyle->PointSize[0.015],
      Points->False, 
      AspectRatio->1/2,
      DefaultFont->{"Times-Roman",14},
      FrameTicks->{None,Automatic},
      Interpolate->True, 
      Causal->False, 
      Stretch->2.0/length,
      InterpStyle->Thickness[0.0], 
      Frame->True,
      options
   ]
]



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


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

(* End[] *) (* SCMTheory`Private` *)


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

(* 
Protect[

]
 *)


(* EndPackage[] *)

(* :Examples:
 *)



