(*From the Continuous Time Fourier Transform to the
Discrete Time Fourier Transform by Sampling
by Nasser M. Abbasi
version April 7 2010*)
Manipulate[
process[Fs, {f1, f2}, {a1, a2}, T, t0, showSinc == 2,
showSamples == 2, unitsInHz == 1, maxCTFTscale, maxDTFTscale,
maxXscale, interpolant, splineOrder],
Panel[Labeled[Grid[{
{Text@Style[Column[{
Row[{Subscript[Style["a", Italic], 1] , "cos(", 2 \[Pi],
Subscript[Style["f", Italic], 1], "(",
Style["t", Italic] - Subscript[Style["t", Italic], 0],
"))"}],
Row[{" + ", Subscript[Style["a", Italic], 2] , "cos(",
2 \[Pi], Subscript[Style["f", Italic], 2], "(",
Style["t", Italic] - Subscript[Style["t", Italic], 0],
"))"}]
}], 10]},
{Control[{{a1, N[9/10],
" \!\(\*SubscriptBox[\(a\), \(1\)]\)"}, -1., 1, 1/10,
ImageSize -> Small, Appearance -> "Labeled" }]},
{Control[{{a2, N[8/10],
" \!\(\*SubscriptBox[\(a\), \(2\)]\)"}, -1., 1, 1/10,
ImageSize -> Small, Appearance -> "Labeled" }]},
{Control[{{f1, N[9/10], " \!\(\*SubscriptBox[\(f\), \(1\)]\)"},
0., 1, 1/10, ImageSize -> Small, Appearance -> "Labeled" }]},
{Control[{{f2, N[2/10], " \!\(\*SubscriptBox[\(f\), \(2\)]\)"},
0., 1, 1/10, ImageSize -> Small, Appearance -> "Labeled" }]},
{Control[{{t0, 0., " \!\(\*SubscriptBox[\(t\), \(0\)]\)"}, 0.,
10, 1/10, ImageSize -> Small, Appearance -> "Labeled" }]}
}, Spacings -> {0, 0}],
Text@Style[
Row[{"signal ", Subscript[Style["x", Italic],
Style["a", Italic]], "(", Style["t", Italic],
") specification"}], 10], {{Top, Center}}], FrameMargins -> 1,
ImageSize -> 190],
Panel[Labeled[
Control[{{T, 6, " T"}, 1, 6, 1, ImageSize -> Small,
Appearance -> "Labeled" }],
Style["signal duration (sec)", 10], {{Top, Center}},
Spacings -> {0, 0}], FrameMargins -> 3, ImageSize -> 190],
Panel[Labeled[
Control[{{Fs, 4.5, " \!\(\*SubscriptBox[\(f\), \(s\)]\)"}, .1, 10,
1/10, ImageSize -> Small, Appearance -> "Labeled"}],
Style["sampling frequency (Hz)", 10], {{Top, Center}},
Spacings -> {0, 0}], FrameMargins -> 3, ImageSize -> 190],
Panel[Labeled[Grid[{
{
Control[{{interpolant, 2, ""}, {1 -> "sinc", 2 -> "staircase",
3 -> "spline"}, ControlType -> RadioButtonBar,
Appearance -> "Vertical"}]
,
Control[{ {splineOrder, 0, Style["order", 10]},
{0 -> Style["0", 9],
1 -> Style["1", 9],
2 -> Style["2", 9],
3 -> Style["3", 9],
4 -> Style["4", 9]
},
Enabled -> Dynamic@TrueQ[interpolant == 3],
ControlType -> PopupMenu, ImageSize -> Tiny}]
}
}, Spacings -> {1, 1}, BaselinePosition -> Bottom,
Alignment -> Bottom],
Style["select algorithm for interpolant", 10], {{Top, Center}}],
Alignment -> {Center, Center}, FrameMargins -> 1, ImageSize -> 190],
Panel[Labeled[Grid[{
{Panel[Labeled[
Control[{{unitsInHz, 2, ""}, {1 -> "Hz", 2 -> "radial"},
ControlType -> RadioButtonBar, Appearance -> "Vertical"}],
Style["x axis units", 10], {{Top, Center}}],
Alignment -> {Center, Center}, FrameMargins -> 1,
ImageSize -> {70, 70}],
Panel[Grid[{{
Labeled[Control[{{showSinc, 2, ""}, {1, 2},
ControlType -> Checkbox}],
Style["show interpolant", 10], {{Top, Center}}],
SpanFromLeft},
{Labeled[
Control[{{showSamples, 1, ""}, {1, 2},
ControlType -> Checkbox}],
Style["show samples", 10], {{Top, Center}}], SpanFromLeft}
}, Spacings -> {0, 0}], Alignment -> {Center, Center},
FrameMargins -> 1, ImageSize -> 110], SpanFromLeft},
{Panel[Labeled[Column[{
Control[{{maxXscale, 7,
" \!\(\*SubscriptBox[\(x\), \(a\)]\)(t)"}, .1, 7, 1/10,
ImageSize -> Tiny, Appearance -> "Labeled" }],
Control[{{maxCTFTscale, 2, " CTFT"}, .1, 3, 1/10,
ImageSize -> Tiny, Appearance -> "Labeled" }],
Control[{{maxDTFTscale, .5, " DTFT"}, .5, 3, 1/10,
ImageSize -> Tiny, Appearance -> "Labeled" }]},
Alignment -> Center
], Style["adjust max x axis scale", 10], {{Top, Center}}
], Alignment -> {Center, Center}, FrameMargins -> 0,
ImageMargins -> 0, ImageSize -> 180], SpanFromLeft}},
Spacings -> {0, 0}
],(*Grid[...]*)
Style["plot options", 9], {{Top, Center}},
Spacings -> {0, 0}],(*Labeled*)
Alignment -> {Center, Center}, FrameMargins -> 0, ImageMargins -> 0,
ImageSize -> 190],(*Panel*)
ControlPlacement -> Left,
FrameMargins -> 0,
ImageMargins -> 0,
ContinuousAction -> False,
SynchronousUpdating -> False,
AutorunSequencing -> {1, 2, 3},
Initialization :>
{
(* The following 3 functions are used by plotting to control ticks*)
xticksForDTFT[unitsInHz_, tickLen_][min_, max_] :=
Module[{},
If[unitsInHz, {{-4, -4, tickLen, Red}, {-3, -3, tickLen,
Red}, {-1.5, -1.5, tickLen, Red}, {-.5, -.5, tickLen,
Red}, {.5, .5, tickLen, Red}, {1.5, 1.5, tickLen, Red}, {3, 3,
tickLen, Red}, {4, 4, tickLen, Red}}, {{-3, -6 Pi, tickLen,
Red}, {-1.5, -3 Pi, tickLen, Red}, {-1/2, -Pi, tickLen,
Red}, {1/2, Pi, tickLen, Red}, {1.5, 3 Pi, tickLen, Red}, {3,
6 Pi, tickLen, Red}}
]
];
ticksForCTFTinRadian[tickLen_][min_, max_] := Module[{},
Table[{i, i*2 Pi, .05, Red}, {i, Ceiling[min], Floor[max], 1}]];
ticksForCTFTinHz[tickLen_][min_, max_] := Module[{},
Table[{i, i, tickLen, Red}, {i, Ceiling[min], Floor[max], 1}]];
(* This is the sinc interpolation formula *)
(* T: Sampling period *)
(* samples: list of samples, sequence of numbers *)
(* t: the instance of time to evaluate interpolant at *)
sincInterpolate[t_, T_, samples_] := Module[{k},
Sum[samples[[ k ]]*Sinc[Pi/T (t - (k - 1) T)], {k, 1,
Length[samples]}]
];
generateSignalExpression[fUser_, a_, t0_] :=
Module[{m, s, xa, xa2, xa1, c},
m = Abs[a[[1]]];
s = Text@If[fUser[[1]] == 0, ToString[m],
Row[{m, " cos(2 \[Pi] \[Times] ", ToString[fUser[[1]]], " (",
Style["t", Italic], " - ", ToString[t0], "))"}]];
If[m != 0, If[a[[1]] < 0, xa1 = Text@Row[{"- ", s}], xa1 = s]];
m = Abs[a[[2]]];
s = Text@If[fUser[[2]] == 0, ToString[m],
Row[{m, " cos(2 \[Pi] \[Times] ", ToString[fUser[[2]]], " (",
Style["t", Italic], " - ", ToString[t0] , "))"}]];
If[m != 0,
If[a[[2]] < 0, xa2 = Text@Row[{" - ", s}],
xa2 = Text@Row[{" + ", s}]]];
Which[
a[[1]] == 0 && a[[2]] == 0, xa = Text["0"],
a[[1]] == 0, xa = xa2,
a[[2]] == 0, xa = xa1,
True, xa = Text@Row[{xa1, xa2}]
];
StringReplace[ToString[xa, FormatType -> TraditionalForm],
c : LetterCharacter ~~ "$" ~~ DigitCharacter .. :> c]
];
(* Main Manipulate expression *)
process[Fs_, fUser_, a_, duration_, t0_, showSinc_, showSamples_,
unitsInHz_, maxCTFTscale_, maxDTFTscale_, maxXscale_,
interpolant_, splineOrder_] :=
Module[{ nSamples, data, xa, ctft, dtft, nHarmonics = 2, xaString,
sincNumberOfPoints = 10, sincDeltaTime, xaReconstruct, pSinc,
panel, disks, lines, dtftTerm, tickLen = .05, from, to,
ticksFontSize = 9, aspectRatio = .49, imageSize = 186,
spectraPlotOptions, k, f, t, i, n, Ts = 1/Fs, p1, p2, p3, p5,
p6, x},
xaString = generateSignalExpression[fUser, a, t0];
xa =
Sum[a[[k]] Cos[
2 Pi fUser[[k]] (t - t0)] (UnitStep[t - t0] -
UnitStep[t - t0 - duration]), {k, 1, nHarmonics}];
nSamples = Floor[duration/Ts] + 1;
(*for speed, I enter the DTFT term directly*)
dtftTerm =
N[a[[1]] Cos[2 fUser[[1]] \[Pi] (n Ts)] +
a[[2]] Cos[2 fUser[[2]] \[Pi] (n Ts)]];
dtft =
Sum[dtftTerm*Exp[-I 2 Pi f n], {n, 0, nSamples}]*
Exp[-I 2 Pi f Floor[t0/Ts]];
(*for speed, I enter the CTFT term directly*)
ctft =
Chop[1/(4 \[Pi])
I (a[[1]] (E^(-2 I \[Pi] (-fUser[[1]] duration +
f (duration + t0)))/(f - fUser[[1]]) +
E^(-2 I \[Pi] (fUser[[1]] duration +
f (duration + t0)))/(f + fUser[[1]]) - (
2 E^(-2 I f \[Pi] t0) f)/(f^2 - fUser[[1]]^2)) +
a[[2]] (E^(-2 I \[Pi] (-fUser[[2]] duration +
f (duration + t0)))/(f - fUser[[2]]) +
E^(-2 I \[Pi] (fUser[[2]] duration +
f (duration + t0)))/(f + fUser[[2]]) - (
2 E^(-2 I f \[Pi] t0) f)/(f^2 - fUser[[2]]^2)))];
from = t0;
to = t0 + maxXscale;
If[to > t0 + duration, to = t0 + duration];
data =
Table[{t0 + (i - 1)*Ts, dtftTerm /. n -> (i - 1)}, {i, 1,
nSamples}];
(* This puts little dots the sample location *)
disks =
Table[{PointSize[.014], Red, Opacity[1],
Point[{data[[i, 1]], data[[i, 2]]}]}, {i, 1, nSamples}];
(* This draws the samples lines *)
lines =
Table[Line[{{data[[i, 1]], 0}, {data[[i, 1]],
data[[i, 2]]}}], {i, 1, nSamples}];
p1 = Plot[xa, {t, from, to}, AxesOrigin -> Automatic,
Ticks -> {Automatic, Automatic},
ImagePadding -> {{24, 45}, {5, 5}},
PlotRange -> {{from, to}, {-2, 2}},
ImageMargins -> 0,
ImageSize -> 380,
AspectRatio -> .40,
AxesLabel -> {Style[Row[{Style["t", Italic], " (sec)"}], 10],
Style[Row[{Subscript[Style["x", Italic], Style["a", Italic]],
"(", Style["t", Italic], ")"}], 12]},
PlotLabel ->
If[Max[fUser]*2 > Fs,
Style["Warning: signal is undersampled!", Red]],
TicksStyle -> Directive[ticksFontSize],
Epilog -> {If[showSamples, {{Dashed, Red, disks}, lines},
Text[" "]]}
];
spectraPlotOptions = {PlotRange -> All,
TicksStyle -> Directive[ticksFontSize], ImageSize -> imageSize,
AxesOrigin -> {0, 0},
AspectRatio -> aspectRatio, Exclusions -> None};
p2 =
Plot[ComplexExpand[Abs[ctft]], {f, -maxCTFTscale, maxCTFTscale},
Evaluate[spectraPlotOptions],
ImagePadding -> {{10, 55}, {25, 20}},
Ticks -> {If[unitsInHz, ticksForCTFTinHz[tickLen],
ticksForCTFTinRadian[tickLen]], Automatic},
AxesLabel -> {Row[{If[unitsInHz, Style["f", 10, Italic],
Style["\[CapitalOmega]", 10]],
Style[If[unitsInHz, " (Hz)", " (rad/sec)"], 9]}],
If[unitsInHz,
Style[Row[{"|", Subscript[Style["X", Italic],
Style["a", Italic]], "(", Style["f", Italic], ")|"}],
12], Style[
"|\!\(\*SubscriptBox[\(X\), \(a\)]\)(\[CapitalOmega])|",
12]]}
];
p3 =
Plot[ComplexExpand[Arg[ctft]], {f, -maxCTFTscale, maxCTFTscale},
Evaluate[spectraPlotOptions],
ImagePadding -> {{10, 55}, {5, 20}},
Ticks -> {If[unitsInHz, ticksForCTFTinHz[tickLen],
ticksForCTFTinRadian[tickLen]], {-Pi, 0, Pi}},
AxesLabel -> {Row[{If[unitsInHz, Style["f", Italic, 10],
Style["\[CapitalOmega]", 10]],
Style[If[unitsInHz, " (Hz)", " (rad/sec)"], 9]}],
If[unitsInHz,
Style[Row[{"phase", Subscript[Style["X", Italic],
Style["a", Italic]], "(", Style["f", Italic], ")"}], 12],
Style[Row[{"phase ", Subscript[Style["X", Italic],
Style["a", Italic]], "(\[CapitalOmega])"}], 12]]}
];
Which[interpolant == 1 || interpolant == 2,
{
Which[interpolant == 1,(*sinc*)
{
sincDeltaTime = Ts/(sincNumberOfPoints - 1);
xaReconstruct = Table[{t + t0,
sincInterpolate[t, Ts, data[[All, 2]] ]}, {t, 0, duration,
sincDeltaTime}
];
},
interpolant == 2,(*stairacse*)
xaReconstruct = Table[{
{(t0 + (i - 1)*Ts) - Ts/2, data[[i, 2]]},
{(t0 + (i - 1)*Ts) + Ts/2, data[[i, 2]]},
{(t0 + (i - 1)*Ts) + Ts/2,
If[i == nSamples, 0, data[[i + 1, 2]]]}},
{i, 1, nSamples}];
];
pSinc = ListPlot[xaReconstruct,
Joined -> True, PlotRange -> All, PlotStyle -> {Magenta},
ImageMargins -> 5, ImageSize -> imageSize];
},
True,(*spline*)
{
If[nSamples > 1,(*can't do spline on one sample*)
pSinc =
Plot[Interpolation[data, InterpolationOrder -> splineOrder][
x], {x, data[[1, 1]], data[[-1, 1]]}, ImageMargins -> 5,
ImageSize -> imageSize, PlotStyle -> {Magenta}],
pSinc = ListPlot[{0, 0}]];
}
];
p5 =
Plot[Abs[dtft], {f, -maxDTFTscale, maxDTFTscale},
Evaluate[spectraPlotOptions],
ImagePadding -> {{10, 50}, {25, 20}},
Ticks -> {xticksForDTFT[unitsInHz, tickLen], Automatic},
AxesLabel -> {Row[{If[unitsInHz, Style["f", 10, Italic],
Style["\[Omega]", 10]],
Style[If[unitsInHz, " (Hz)", " (rad)"], 9]}],
If[unitsInHz,
Style[Row[{"|", Style["X", Italic], "(", Style["f", Italic],
")|"}], 12],
Style[Row[{"|", Style["X", Italic], "(\[Omega])|"}], 12]]}
];
p6 =
Plot[Arg[dtft], {f, -maxDTFTscale, maxDTFTscale},
Evaluate[spectraPlotOptions],
Ticks -> {xticksForDTFT[unitsInHz, tickLen], {-Pi, Pi}},
ImagePadding -> {{10, 50}, {5, 20}},
AxesLabel -> {Row[{If[unitsInHz, Style["f", Italic, 10],
Style["\[Omega]", 10]],
Style[If[unitsInHz, " (Hz)", " (rad)"], 9]}],
If[unitsInHz,
Style[Row[{"phase ", Style["X", Italic], "(",
Style["f", Italic], ")"}], 12],
Style[Row[{"phase ", Style["X", Italic], "(",
Style["\[Omega]", Italic], ")"}], 12]]}
];
panel = Panel[Grid[{
{Null, SpanFromLeft},
{Text@
Row[{Subscript[Style["x", Italic, 12], Style["a", Italic]],
"(", Style["t", Italic, 12], ") = ",
Style[xaString, 12]}], SpanFromLeft},
{Null,
SpanFromLeft}, {Text@
Style[Row[{"Nyquist frequency = ",
NumberForm[Max[fUser]*2, {3, 2}], " Hz"}], 12],
Text@Style[
Row[{"sampling period ", Subscript[Style["T", Italic],
Style["s", Italic]], " = ", NumberForm[Ts, {3, 2}],
" sec"}], 12]}
}, Frame -> None, Alignment -> Center, Spacings -> {6, 0}],
FrameMargins -> 1, ImageMargins -> 0];
Grid[
{
{ Labeled[ If[showSinc, Show[{p1, pSinc}], p1],
panel, {{Top, Center}}, Spacings -> {0, 1}], SpanFromLeft},
{Labeled[Grid[{{p2, p3}}],
Text@Style["CTFT, continuous\[Hyphen]time Fourier transform",
12, Bold], {{Top, Center}}, Spacings -> {0, 0}],
SpanFromLeft},
{Labeled[Grid[{{p5, p6}}],
Text@Style["DTFT, discrete\[Hyphen]time Fourier transform",
12, Bold], {{Top, Center}}, Spacings -> {0, 0}],
SpanFromLeft}
},
Alignment -> Center,
Spacings -> {0, 1},
Frame -> All
]]}]