(* :Version: Mathematica 2.0 *)

(* :Name: FMPlot.m *)

(* :Title: Frequency Modulation Spectral Plotting *)

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

(* :Keywords: Frequency Modulation, FM Synthesis, Sound, Music *)

(* :Requirements: none. *)

(* :Warnings: none. *)

(* :Date: 22 March 1995 *)

(* :Summary:

This package contains functions for plotting and manipulating
the frequency spectra of various carrier/frequency/index combinations
used in simple FM synthesis.  In addition, additive-synthesis 
spectra can be created and manipulated with this package.  There
is a tutorial that goes along with this package called FMPlot-tutorial.ma.

*)

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



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


AbsPositiveFMSpectrum::usage = "AbsPositiveSpectrum[carrierFrequency, \
   modulatorFrequency, index, minimumAmplitude:0.045] returns the \
   positive frequencies created by the particular FM combination with \
   an absolutes value for the amplitudes.";

CombineSpectra::usage = "CombineSpectra[Spectrum1, Spectrum2, Spectrum3, ...] \
   adds the given spectra into one spectal list which is returned.";

Deviation::usage = "Deviation[modulatorFrequency, index] returns the \
   deviation given the modulatorFrequency and index.";
  
FMPlot::usage = "FMPlot[carrierFrequency, modulatorFrequency, \
   index, minimumAmplitude:0.045, width:0.02, plotOptions] plots a \
   frequency spectra of the given carrier / modulator / index; where \
   width controls the size of the spectral lines, and lowering the \
   minimumAmplitude may extend the range of the plot to include less \
   important sidebands.";

FMSpectrum::usage = "FMSpectrum[carrierFrequency, modulatorFrequency, index, \
   minimumAmplitude:0.045] returns the partials created by the particular FM \
   combination.";

Index::usage = "Index[modulatorFrequency, deviation] returns the \
   index given the modulatorFrequency and deviation.";

NormalizeSpectrum::usage = "NormalizeSpectrum[list, max:1] normalizes the \
   second number in a list of number pairs to the value of max which is \
   defaulted to 1.0";
  
PlotAbsPositiveFMSpectrum::usage =
   "PlotAbsPositiveFMSpectrum[carrierFrequency, modulatorFrequency, \
   index, minimumAmplitude:0.045, width:0.02, plotOptions] plots a \
   frequency spectra of the given carrier / modulator / index; where \
   width controls the size of the spectral lines, and lowering the \
   minimumAmplitude may extend the range of the plot to include less \
   important sidebands.";

PositiveFMSpectrum::usage = "PositiveFMSpectrum[carrierFrequency, \
   modulatorFrequency, index, minimumAmplitude:0.045] returns the positive \
   frequencies created by the particular FM combination.";

RemovePhase::usage = "RemovePhase[spectralList] removes phase information \
   from the spectralList.";

SpectralAmp::usage = "SpectralAmp[spectralList] returns a list of \
   only the amplitudes in the spectral list.";

SpectralFreq::usage = "SpectralFreq[spectralList] returns a list of \
   only the frequencies in the spectral list.";

SpectralLines::usage = "SpectralLines[spectralList,  width:0.02] gives the \
   graphs for a set of spectral lines of width thickness.";

SpectralListQ::usage = "SpectralListQ[spectralList] tests whether the \
   argument is a well-formed spectral list.";

SpectralPlot::usage = "SpectralPlot[spectralList, lineWidth_:0.02] plots a \
   frequency spectra of the given list of spectral lines of the form \
   {{frequency ,amplitude},...}.";



Begin["`Private`"]



(***** AbsPositiveFMSpectrum *********************************************)

AbsPositiveFMSpectrum[
   carrier_?NumberQ,  modulator_?NumberQ,
   index_ /; index >= 0,  min_:0.045] :=
      Map[Abs, PositiveFMSpectrum[carrier,modulator,index,min]]



(***** CombineSpectra ****************************************************)

CombineSpectra[spectra___] :=
   Module[
      {
         list,
         tempList,
         tempLength,
         i,j
      },
 
      tempList = Sort[Flatten[{spectra}, 1]];

      (* create the bins for frequency amplitude info *)
      
      list = Apply[Union, Map[List, Transpose[
         {1,0} Transpose[tempList]]]];

      (* fill the frequency bins with amplitude info *)

      tempLength = Length[tempList];
      list[[1,2]] += tempList[[1,2]];
      For[i=2; j=1, i<=tempLength, i++,
         If[tempList[[i,1]] != tempList[[i-1,1]], j++];
         list[[j,2]] += tempList[[i,2]]
      ];

      list
   ]



(***** Deviation *********************************************************)

Deviation[modulator_?NumberQ, index_ /; index >= 0] := modulator index



(***** FMPlot ************************************************************)

(* FMPlot is the same function as PlotAbsPositiveFMSpectrum *)

FMPlot[
   carrier_?NumberQ, modulator_?NumberQ, 
   index_ /; index >= 0, min_:0.045, width_:0.02, opts___] :=
      PlotAbsPositiveFMSpectrum[carrier, modulator, index, min, width, opts]



(***** FMSpectrum ********************************************************)

FMSpectrum[
   carrier_?NumberQ, modulator_?NumberQ, 
   index_ /; index >= 0, min_:0.045] :=
   Module[
      {
         k = 1,
         i = index,
         list
      },
      
      list = {{carrier, J[0,i]}};
      
      While[ (J[k,i] > min) || (k <= Round[index]+2) ,
         AppendTo[list,{carrier + k modulator, J[k,i]}];
         AppendTo[list,{carrier - k modulator, J[k,i]}];
         k++;
      ];
      Sort[list]
   ]



(***** Index *************************************************************)

Index[modulator_?NumberQ, deviation_?NumberQ] := deviation / modulator


      
(***** NormalizeSpectrum *************************************************)

NormalizeSpectrum[spectralList_?SpectralListQ, goal_:1.0] :=
   Module[ {scaler = N[goal/Max[SpectralAmp[spectralList]]]},
      Transpose[ {1, scaler} Transpose[spectralList] ]
   ]



(***** PlotAbsPositiveFMSpectrum *****************************************)

PlotAbsPositiveFMSpectrum[
   carrier_?NumberQ, modulator_?NumberQ, 
   index_ /; index >= 0, minimumAmplitude_:0.045, width_:0.02, opts___] :=
   Module[
      {
         titleForm,
         title = "Carrier: "     <> ToString[carrier]   <>
                 ", Modulator: " <> ToString[modulator] <>
                 ", Index: "     <> ToString[index],
         spectrum = NormalizeSpectrum[AbsPositiveFMSpectrum[carrier,
            modulator, index, minimumAmplitude]],
      },
      titleForm = FontForm[title,{"Palatino-Italic",16}];

      Show[
         ListPlot[spectrum,
            opts,
            PlotStyle->{Thickness[0.0],RGBColor[1,1,1]},
            PlotJoined->False,
            Frame->True,
            AxesOrigin->{0,0},
            DefaultFont->{"Helvetica-Bold", 12},
            FrameLabel->{"Frequency [Hz]", "Amplitude [normalized]",
               titleForm, " "},
            DisplayFunction->Identity,
            PlotRange->All
         ],
         SpectralLines[spectrum, width],
         DisplayFunction :> $DisplayFunction,
         PlotRange->All
      ]
   ]



(***** PositiveFMSpectrum ************************************************)

PositiveFMSpectrum[
   carrier_?NumberQ, modulator_?NumberQ, 
   index_ /; index >= 0, min_:0.045] :=
   Module[{list},
      list = FMSpectrum[carrier,modulator,index,min];
      list = Sort[Map[If[#[[1]]<0, -#, #]&, list]];
      CombineSpectra[list]  
   ]



(***** RemovePhase *******************************************************)

RemovePhase[spectralList_?SpectralListQ] :=
   Module[{list},
      list = Sort[Map[If[#[[1]]<0, -#, #]&,spectralList]];
      Map[(Abs[#])&,CombineSpectra[list]]
   ]



(***** SpectralAmp *******************************************************)

SpectralAmp[spectralList_?SpectralListQ] :=
   Transpose[spectralList][[2]]



(***** SpectralFreq ******************************************************)

SpectralFreq[spectralList_?SpectralListQ] :=
   Transpose[spectralList][[1]]



(***** SpectralLines *****************************************************)

SpectralLines[spectrum_?SpectralListQ, width_:0.02] :=
   Module[
      {
         axisList,
         lineSet
      },
      axisList = Transpose[{1,0} Transpose[spectrum]];
      lineSet = Transpose[{axisList, spectrum}];
      Graphics[{Thickness[width],Map[Line,lineSet]}]
   ]   



(***** SpectralListQ *****************************************************)

SpectralListQ[spectralList_] :=
   If[ListQ[spectralList],
      Apply[And, Map[(Length[#] == 2)&,spectralList] ],
      False
   ]



(***** SpectralPlot ******************************************************)

SpectralPlot[spectralList_?SpectralListQ, width_:0.02, opts___] :=
   Module[{extractedList = spectralList},

      (* special features so that animation can be done - see tutorial *)

      If[(First[Rest[extractedList]][[2]] === 0) &&
         (First[extractedList][[2]] != 0),
         extractedList = Drop[extractedList,1]
      ];
      If[First[extractedList][[2]] === 0,
         extractedList = Drop[extractedList,1]
      ];
      If[(extractedList[[Length[extractedList]-1]][[2]] === 0) &&
         (Last[extractedList][[2]] != 0),
         extractedList = Drop[extractedList,-1]
      ];
      If[Last[extractedList][[2]] === 0,
         extractedList = Drop[extractedList,-1] 
      ];

      Show[
         ListPlot[spectralList,
            opts,
            PlotStyle->{Thickness[0.0],RGBColor[1,1,1]},
            PlotJoined->False,
            Frame->True,
            AxesOrigin->{-100000,0},
            DefaultFont->{"Helvetica-Bold", 12},
            FrameLabel->{"Frequency", "Amplitude", 
               FontForm["'",{"Courier", 2}], " "},
            DisplayFunction->Identity,
            PlotRange->All
         ],
         SpectralLines[extractedList, width],
         DisplayFunction :> $DisplayFunction,
         PlotRange->All
      ]
   ]



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

J[k_, n_]:= N[ BesselJ[k,n] ]



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

End[ ]   (* FMPlot`Private` *)


Protect[

   AbsPositiveFMSpectrum,
   CombineSpectra, 
   Deviation, 
   FMPlot,
   FMSpectrum, 
   Index, 
   NormalizeSpectrum, 
   PlotAbsPositiveFMSpectrum,
   PositiveFMSpectrum, 
   RemovePhase,
   SpectralAmp,
   SpectralFreq,
   SpectralLines, 
   SpectralListQ, 
   SpectralPlot,

]

EndPackage[]


(* :Examples:
 
   See also the tutorial.

   Spectrum from a carrier=1, modulator=2, and index=4:
      FMSpectrum[1,2,4]
   
   Now fliping all negative frequencies into the positive domain:
      PositiveFMSpectrum[1,2,4]
   
   Now removing phase information from the spectrum:
      Spectrum124 = AbsPositiveFMSpectrum[1,2,4]
   
   Now  normalizing the spectrum amplitudes:
      NormalizeSpectrum[spectrum124]
   
   Or normalize the spectrum as percentages:
      Normalize[spectrum124,100]
   
   Example Plot:
      FMPlot[200,50,1]
   
   Plot with lowered threshold of amplitude:
      FMPlot[200, 50, 1, 0.000001]
   
   Or with a different thickness to the spectral lines:
      FMPlot[200, 50, 1, 0.000001, 0.01]
   
      FMPlot[440,880,20]
      FMPlot[1000,100,40,0.045,0.005]

      SpectralLines[FMSpectrum[1,2,3]] // Show
   
 *)

