(*Power content of angle modulation (FM and PM)
by Nasser M. Abbasi
version September 6 2009*)

Manipulate[
 If[typeOfModulation == FM, k = Max[Min[k, 2], 0], 
  k = Max[Min[k, 3.5], 0]];
 
 frequencyModulation[k, fc, fm, typeOfModulation, Am, Ac],
 
 Grid[{
   {
    Item[Control[{{typeOfModulation, FM, 
        "select type of modulation"}, {FM -> "FM", PM -> "PM"}, 
       ControlType -> SetterBar, ImageSize -> Tiny}]
     , Alignment -> Left],
    
    Item[Control[{{k, 1.14, "k, deviation constant "}, 0, 
       Dynamic[If[typeOfModulation == FM, 2, 3.5]], .01, 
       Appearance -> "Labeled", ImageSize -> Tiny}]
     , Alignment -> Right]
    },
   
   {
    Item[Control[{{Am, 1, 
        "\!\(\*SubscriptBox[\(A\), \(m\)]\), message amplitude "}, 
       0.5, 1.5, .01, Appearance -> "Labeled", ImageSize -> Tiny}]
     , Alignment -> Right],
    
     Item[
     Control[{{fc, 6.7, 
        "\!\(\*SubscriptBox[\(f\), \(\(c\)\(,\)\(\\\ \)\)]\)carrier \
frequency (hz)"}, 5, maxCarrierFreq, .01, Appearance -> "Labeled", 
       ImageSize -> Tiny}]
     , Alignment -> Right]
    },
   
   {
    Item[Control[{{Ac, 1, 
        "\!\(\*SubscriptBox[\(A\), \(c\)]\), carrier amplitude "}, .5,
        1.5, .01, ImageSize -> Tiny, Appearance -> "Labeled"}]
     , Alignment -> Right],
    
    Item[Control[ {{fm, 0.88, 
        "\!\(\*SubscriptBox[\(f\), \(\(m\)\(,\)\(\\\ \)\)]\)message \
frequency (hz)"}, .18, maxMessageFreq, .01, Appearance -> "Labeled", 
       ImageSize -> Tiny}]
     , Alignment -> Right]
    }
   },
  Spacings -> {2, 0},
  Frame -> None
  ],
 
 ControlPlacement -> Top,
 SynchronousUpdating -> False,
 AutorunSequencing -> Automatic,
 Initialization :>
  (
   maxCarrierFreq = 8;
   maxMessageFreq = 1;
   tmax = 4;
   fontSizeForTitles = 14;
   fontSizeForSubTitles = 11;
   FM = 1;
   PM = -1;
   
   (********************************************************)
   (* This is the main function called by Manipulate[]     *)
   (********************************************************)
   frequencyModulation[k_, fc_, fm_, typeOfModulation_, Am_, Ac_] := 
    Module[{message, carrier, wm, wc, msg, \[Beta]},
     
     wm = 2 Pi fm;
     wc = 2 Pi fc;
     carrier = Ac Cos[wc t];
     
     (*set \[Beta], the modulation index, 
     depending on the modulation type*)
     If[typeOfModulation == 
       FM, {message = Am Cos[wm t], \[Beta] = k Am/fm}, {message = 
        Am Sin[wm t], \[Beta] = k Am}];
     
     GraphicsGrid[{
       {
        msg =
         Style[
          "modulation index \[Beta] = " <> 
           ToString[ NumberForm[\[Beta] , {3, 2}]], 
          fontSizeForSubTitles];
        
        plotLabel = 
         Style[Column[{ 
            Style[If[typeOfModulation == FM, "FM modulated carrier", 
              "PM modulated carrier"], fontSizeForTitles] ,
            msg}
           , Center]] ;
        
        Plot[{Ac Cos[wc t + \[Beta]   Sin[wm t] ], message}, {t, 0, 
          tmax },
         Ticks -> {Automatic, Automatic},
         AxesLabel -> {Style[Row[{Style["t", Italic], " (sec)"}], 
            Larger]}, PlotLabel -> plotLabel,
         PlotStyle -> {Red, Blue},
         AspectRatio -> 1/4
         ],
        SpanFromLeft
        },
       {
        frequencyPlot[fm, fc, \[Beta], Ac, Am], SpanFromLeft
        }
       ,
       {
        powerRatioPlot[fm, \[Beta]], SpanFromLeft
        }
       },
      Spacings -> {Scaled[.1], Scaled[.1]},
      Frame -> {False, All},
      Dividers -> None,
      AspectRatio -> Full,
      ImageSize -> {520, 390}]
     ];
   
   (********************************************************)
   (* This function generates the plot for the power ratio *)
   (********************************************************)
   powerRatioPlot[fm_, \[Beta]_] := 
    Module[{bandwidth, nSideTermsToUse, data, bessel, tmp, n, xticks, 
      p1, p2, colors, plotLabel},
     
     bandwidth = 2 (\[Beta] + 1) fm;
     nSideTermsToUse = Floor[bandwidth/(2*fm)] + 1;
     tmp = BesselJ[0, \[Beta]]^2;
     bessel = Table[BesselJ[n, \[Beta]]^2, {n, 1, nSideTermsToUse}];
     data = 
      Table[{n, tmp + 2 Total[bessel[[1 ;; n]]]}, {n, 1, 
        nSideTermsToUse}];
     data = Insert[data, {0, tmp}, 1];
     
     xticks = Table[{n, ToString[2 n]}, {n, 0, nSideTermsToUse}];
     
     (*I do not know how to make lines joined and also make filling*)
\

     (*as vertical lines using the same options as they conflict    *)
\

     (*so I make 2 plots one with vertical lines, 
     and one joined   *)
     
     colors = 
      Table[If[
        n == 0, {Red, PointSize[.015]}, {Blue, PointSize[.015]}], {n, 
        0, nSideTermsToUse}];
     
     plotLabel = 
      Style[Column[{ 
         Style["power content (normalized) as a function of \
bandwidth", fontSizeForTitles] },
        , Center]];
     
     p1 = Show[ListPlot[{#[[2]]},
          PlotLabel -> plotLabel ,
          AspectRatio -> 1/4,
          Ticks -> {xticks, {0, 0.25, 0.5, 0.75, 1}},
          AxesOrigin -> {0, 0},
          Filling -> Axis,
          PlotStyle -> #[[1]],
          FillingStyle -> #[[1, 1]],
          
          AxesLabel -> {Style[
             Column[{"number of", "sidbands"}, Center], Larger] },
          Frame -> False,
          Joined -> False] & /@ Transpose[{colors, data}],
       PlotRange -> All];
     
     p2 = ListPlot[data,
       Joined -> True,
       PlotStyle -> {Black, Dashed, Thin}];
     
     Show[p1, p2]
     ];
   
   (********************************************************)
   (* This function generates the plot for the spectra     *)
   (********************************************************)
   frequencyPlot[fm_, fc_, \[Beta]_, Ac_, Am_] := 
    Module[{nSideTermsToUse, yValues, xValues, g, plotLabel, gLabels, 
      roundedYValues, bandwidth},
     
     bandwidth = 2 (\[Beta] + 1) fm;
     nSideTermsToUse = Round[bandwidth/(2*fm)] + 1;
     
     yValues = (Ac/2) BesselJ[
        Range[-nSideTermsToUse, nSideTermsToUse], \[Beta]];
     
     (*flip negative phase indicated by some BesselJ results will be \
nagtive*)
     Table[
      If[yValues[[i]] < 0, yValues[[i]] = yValues[[i]]*-1, 
       yValues[[i]]], {i, 1, Length[yValues]}];
     
     xValues = 
      Table[fc + n*fm, {n, -nSideTermsToUse, nSideTermsToUse}];
     
     g = Table[{ If[n == nSideTermsToUse + 1,  Red, Blue]   , 
        Line[{{xValues[[n]], 0}, {xValues[[n]], yValues[[n]]}}]}, {n, 
        1, Length[xValues]}];
     
     roundedYValues = Round[yValues, 0.01];
     
     gLabels = 
      Table[Text[ 
        Style[ NumberForm[roundedYValues[[n]] , {3, 2}], 
         Small], {xValues[[n]], yValues[[n]] + 0.1 yValues[[n]]}, {-1,
          0}, {0, 1}], {n, 1, Length[xValues]}];
     
     plotLabel = 
      Column[{ 
        Style["modulated carrier spectra (magnitude)", 
         fontSizeForTitles] ,
        Row[{Style[
           "bandwidth (98% of power) = " <> 
            ToString[ NumberForm[bandwidth, {3, 2}]] <> " hz", 
           fontSizeForSubTitles],
          "     ",
          
          Style["number of sidbands = " <> 
            ToString[  2 nSideTermsToUse ], fontSizeForSubTitles]}]}
       , Center] ;
     
     Show[ListPlot[{{fc, 0}},
       PlotLabel -> plotLabel,
       AxesOrigin -> {fc, 0},
       Axes -> {True, False},
       PlotRange -> {Automatic, {-.1 Max[yValues], 1.7*Max[yValues]}},
       AspectRatio -> 1/4,
       AxesLabel -> {Style[Row[{Style["f", Italic], " (hz)"}], 
          Larger]},
       Ticks -> {Automatic, None}
       ],
      Graphics[g], Graphics[gLabels]]
     ]
   )
 ]