(* :Version: Mathematica 2.0 *)

(* :Name: TDAC.m *)

(* :Title: Time Domain Alias Cancellation Package *)

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

(* :Keywords: Audio Compression, Time Domain Alias Cancellation *)

(* :Requirements: none. *)

(* :Warnings: none. *)

(* :Date: 7 February 1997 *)

(* :Summary:

This package contains functions to demonstrate the technique
used in audio compression that is called "Time-domain aliasing
cancellation".  See the tutorial for examples.

*)

(* 
 * Changing $Context to TDAC` and $ContextPath to {TDAC`, System`}:
 *)
 
BeginPackage["TDAC`"]



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

IMDCT::usage = "IMDCT[spectrumList] returns the Inverse Modified Discrete \
   Cosine Transform of the spectrum list. The length of spectrumList \
   should be a power of 2 in order to implement the FFT, otherwise the \
   direct DFT will be used to calculate the transform.";

ITDAC::usage = "ITDAC[sequence, frameLength] reverses the process of the \
   TDAC function to give the original signal.  frameLength is the frame \
   size used to generate the sequence with the TDAC function.";

MDCT::usage = "MDCT[signalList] returns the Modified Discrete Cosine \
   Transform of the input list. modSequence should be Real, and the output \
   is real.  Output spectrum length will be 1/2 the length of the input \
   since the full spectrum is symmetric.  The length of the signal list \
   should be a power of 2 in order to implement the FFT, otherwise the \
   direct DFT will be used to calculate the MDCT.";

Resynth::usage = "    ";

SineWindow::usage = "SineWindow[length] returns the first half period \
   of a sampled sinewave shifted by 1/2 of a sample.  Alias is SineW.";

SineW::usage = SineWindow::usage;

TDAC::usage = "TDAC[sequence, frameLength] takes an input sampled signal \
   and returns its modified discrete cosine transform in a form which \
   can then be read by ITDAC to generate the original signal exactly. \
   This function splits the input into frames, zero pads by 1/2 \
   windowLength the beginning and ending of the modSequence. The function \
   adds frameLength/2 extra values to the spectrum output due to the zero \
   padding.frameLength should be a power of 2, and modSequence should be \
   a multiple of frameLength/2 in length. TDAC stands for Time Domain \
   Aliasing Cancellation.";


(* Options *)

Options[Resynth] := {OrigPlot->True, FramePlots->False,
   WindowedFramePlots->False, AnalysisPlot->False,
   ResynthPlots->False, OutputPlot->True}



Begin["`Private`"]



(***** IMDCT *************************************************************)

IMDCT[parameters___] := FastIMDCT[parameters]


FastIMDCT[in_] := Module[
   {n0, z, Z, len = Length[in]*2},
   n0 = N[(len/2+1)/2];
   fullSpectrum = Join[in, -Reverse[in]];
   Z = Chop[N[Table[fullSpectrum[[k+1]]E^(I 2 Pi k n0/len), {k, 0, len-1}]]];
   z = SCMTheory`NumericIDFT[Z];  (* a length of 2^x will use fft *)
   Chop[N[Table[2 (Re[z[[n+1]]] Cos[Pi(n+n0)/len] -
      Im[z[[n+1]]] Sin[Pi (n+n0)/len]), {n,0,len-1}]]]
]


SlowIMDCT[input_] := Module[
   {n0,len = Length[input]*2},
   n0 = N[(len/2+1)/2];
   fullSpectrum = Join[input, -Reverse[input]];
   Chop[N[2/len Table[
      Sum[fullSpectrum[[k+1]] Cos[2Pi/len(n+n0)(k+0.5)],
         {k,0,len-1}], 
      {n,0,len-1}]]]
]


(***** MDCT **************************************************************)

MDCT[parameters___] := FastMDCT[parameters]


FastMDCT[input_] := Module[
   {n0, y, Y, len = Length[input]},
   n0 = N[(len/2+1)/2];
   y = N[Table[input[[n+1]] E^(-I Pi n / len), {n, 0, len-1}]];
   Y = SCMTheory`NumericDFT[y];  (* a length of 2^x will use fft *)
   N[Table[ Re[Y[[k+1]]] *Cos[2Pi n0(k+0.5)/len] + Im[Y[[k+1]]] *
      Sin[2Pi n0(k+0.5)/len], {k, 0, len/2-1}]]
]


DirectMDCT[input_] := Module[ 
   {n0,len = Length[in]},
   n0 = N[(len/2+1)/2];
   N[Table[Sum[in[[n+1]] * Cos[2Pi/len (n+n0)(k+0.5)],
      {n,0,len-1}], {k,0,len-1}]]
]



(***** Resynth ***********************************************************)

Resynth[sig_, wlen_, options___] := Module[{},
   If[check[OrigPlot, options],
      SCMTheory`SeqPlot[sig, Lines->False, Points->False,
         plotStyle["Original Signal"]]];
   frames = Partition[Flatten[Join[{{Zeros[wlen/2]}},
   Map[{#,#}&,Partition[sig, wlen/2]], {{Zeros[wlen/2]}}]], wlen];
   If[check[FramePlots, options], 
      Block[{$DisplayFunction=Identity},
         framePlots = MapIndexed[
            SCMTheory`SeqPlot[#1, Lines->False, Points->False, 
               DefaultFont->{"Times-Roman", 14}, 
               plotStyle["Frame " <> ToString[First[ #2]-1]]]&,
            frames];
      (* framePlots = Join[Partition[framePlots, 2], {{Last[framePlots]}}];*)
      ];
      Show[GraphicsArray[framePlots,GraphicsSpacing->0]]
   ]; 


   If[check[WindowedFramePlots, options],
      Block[{$DisplayFunction=Identity},
         framePlots = MapIndexed[
            SCMTheory`SeqPlot[#1 * SineWindow[wlen], Lines->False, 
            Points->False, DefaultFont->{"Times-Roman", 14}, 
            plotStyle["Windowed Frame " <> ToString[First[#2]-1]]]&, frames];
            SCMTheory`SeqPlot[#1 * SineWindow[wlen], 
               Lines->False, Points->False,
               DefaultFont->{"Times-Roman", 14}, 
        (* framePlots = Join[Partition[framePlots, 2], {{Last[framePlots]}}];*)
      ];
      Show[GraphicsArray[framePlots,GraphicsSpacing->0]]
   ]; 

   analyses = Map[FastMDCT[ SineWindow[wlen]*#]&, frames];
   If[check[AnalysisPlot, options],   
   SCMTheory`SeqPlot[Drop[Drop[TDAC[sig, wlen], -wlen/2], wlen/2],
      AspectRatio->1/4,  Lines->False, Points->False,
      plotStyle["Analysis"]] ];

   reAnalyses = Map[IMDCT, analyses];
   If[True (* check[ResynthPlots], *),
      Block[{$DisplayFunction=Identity},
         framePlots = MapIndexed[
            SCMTheory`SeqPlot[#1 , Lines->False, Points->False, 
            DefaultFont->{"Times-Roman", 14}, 
            plotStyle["Alias Frame " <> ToString[First[#2]-1]]]&, reAnalyses];
         (* framePlots = Join[Partition[framePlots, 2], {{Last[framePlots]}}];*)
      ];
      Show[GraphicsArray[framePlots,GraphicsSpacing->0]]
   ];

   reAnalyses = Map[# * SineWindow[wlen]&, reAnalyses];
   If[True (* check[ResynthPlots], *),
      Block[{$DisplayFunction=Identity},
         framePlots = MapIndexed[
            SCMTheory`SeqPlot[#1 , Lines->False, Points->False,
               DefaultFont->{"Times-Roman", 14}, 
               plotStyle["Windowed Alias " <> ToString[First[#2]-1]]]&,
               reAnalyses];
         (* framePlots = Join[Partition[framePlots, 2], {{Last[framePlots]}}];*)
      ];
      Show[GraphicsArray[framePlots,GraphicsSpacing->0]]
   ];

   subframes =  Drop[ Drop[Partition[Flatten[reAnalyses], wlen/2],1], -1];
   output =  Chop[N[Flatten[ Map[#[[1]]+#[[2]]&, Partition[subframes, 2]]]]];
   If[check[OutputPlot, options], 
      SCMTheory`SeqPlot[output, Lines->False, Points->False, 
         plotStyle["Resynthesized Signal"]]];
   output

*)

]



(***** SineWindow ********************************************************)

SineWindow[wlen_] :=  
   Table[N[Sin[Pi/wlen (n+0.5)]], {n,0, wlen-1}]

SineW[parameters___] := SineWindow[parameters]



(***** ITDAC *************************************************************)

ITDAC[in_, wlen_] := Module[
   {analyses, resynthes, subframes},
   analyses = Partition[in, wlen/2];
   resynthes = Map[# * SineWindow[wlen]&,  Map[IMDCT, analyses] ];
   subframes =  Drop[ Drop[Partition[Flatten[resynthes], wlen/2],1], -1];
   Chop[N[Flatten[ Map[#[[1]]+#[[2]]&, Partition[subframes, 2]]]]]
]



(***** TDAC **************************************************************)

TDAC[in_, wlen_] := Module[
   {frames, analyses},
   frames = Partition[Flatten[Join[{{zeros[wlen/2]}},
   Map[{#,#}&,Partition[in, wlen/2]], {{zeros[wlen/2]}}]], wlen];
   analyses = Map[MDCT[ SineWindow[wlen]*#]&, frames];
   Flatten[analyses] 
]



(***** Private Functions *************************************************)

check[option_, options___] := (option /. {options} /.  Options[Resynth])

plotStyle[title_] := Sequence[Frame->True,
   Lines->False, DefaultFont->{"Times-Roman", 18},
   Points->False,
   PointStyle->PointSize[0.0], 
   FrameLabel->{"","",FontForm[title, {"Times-Bold", 24}],""},	
   FrameTicks->{Automatic, Automatic, None, None}, 
   InterpStyle->Thickness[0.003], 
   Interpolate->True, PadFactor->1]

Zeros[num_] := Table[0, {x,1, num}]



(***** End of Functions **************************************************)

End[ ]   (* TDAC`Private` *)


Protect[

   IMDCT,
   ITDAC,
   MDCT,
   Resynth,
   SineWindow,
   SineW,
   TDAC,

]

EndPackage[]


(* :Examples:
 
   See also the tutorial.

   
 *)

