{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
267 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 
{CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 
0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Time
s" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 26 "Discrete Fourier Transfor
m" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "We d
evelop the idea of a Fourier transform in order to analyze functions i
n creative ways:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 58 "temporal signals in the frequency domain (signal analysis
)" }}{PARA 0 "" 0 "" {TEXT -1 62 "spatial information in the momentum \+
domain (quantum mechanics)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 462 "For some functions the Fourier transform (FT) ca
n be calculated analytically, but we will be mostly concerned with dis
crete FTs for which fast algorithms (FFTs) exist. The DFT (or FFT) is \+
a numerical technique that can be applied to virtually any function. M
any digital measuring devices (oscilloscopes) have built-in hardware w
hich can carry out FFTs for signal analysis. In numerical analysis FFT
s are used to solve ordinary and partial differential equations." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "restart; with(LinearAlgebra)
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "with(plots): with(intt
rans):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "fx:=(4/sqrt(5*Pi)
)/(1+x^2)^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "int(fx^2,x=
-infinity..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "fp
:=fourier(fx,x,p);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "int(f
p^2,p=-infinity..infinity);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 30 "Th
e FT defined as part of the " }{TEXT 19 9 "inttrans " }{TEXT -1 82 "pa
ckage misses a factor, and we redefine the answer so that the norm is \+
preserved." }}{PARA 0 "" 0 "" {TEXT -1 55 "The naive calculation of th
e FT for our function fails:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 59 "fpn:=1/sqrt(2*Pi)*int(fx*exp(I*p*x),x=-infinity..infinity);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "It can be fixed, however, by assum
ing that the variable " }{TEXT 19 1 "p" }{TEXT -1 16 " is real-valued.
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "assume(p,real);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "fpn:=1/sqrt(2*Pi)*int(fx*exp
(I*p*x),x=-infinity..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 57 "fpn:=1/sqrt(2*Pi)*int(fx*cos(p*x),x=-infinity..infinity);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "fpnI:=1/sqrt(2*Pi)*int(fx*si
n(p*x),x=-infinity..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 51 "plot([fp/sqrt(2*Pi),fpn],p=-5..5,color=[red,blue]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 41 "The graph of the expressions is th
e same." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "simplify(fp/sqrt
(2*Pi) - fpn);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "subs(p=3,
%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "simplify(%);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT 257 11 "Exercise 1:" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 26 "Modify the above function
 " }{TEXT 19 2 "fx" }{TEXT -1 186 " in such a way as to change the sca
le over which it varies by replacing the 1 in the denominator by some \+
other constant. Keep the function normalized by adjusting the factor i
n front of " }{TEXT 19 2 "fx" }{TEXT -1 80 ". Calculate the FT, and ob
serve the relationship. Is it true that distributions " }{TEXT 19 2 "f
x" }{TEXT -1 106 " which have more concentration at smaller x-values h
ave transforms that emphasize more the large-p region?" }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "fx:=(2
/Pi)^(1/4)/sqrt(a)*exp(-(x/a)^2);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 16 "fourier(fx,x,p);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 12 "assume(a>0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 16 "fourier(fx,x,p);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "The ass
umption on " }{TEXT 19 1 "a" }{TEXT -1 144 " was required for the FT t
o be meaningful. The example of a Gaussian is particularly illustrativ
e, because it is obvious how the distance scale " }{TEXT 19 1 "a" }
{TEXT -1 49 " appears as the inverse momentum scale in the FT." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "Now let u
s look at examples relevant for signal theory:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 20 "fourier(sin(t),t,w);" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 22 "fourier(cos(3*t),t,w);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 157 "We obtain a Dirac delta at both positive and negative \+
frequency. For signals made up of a few frequencies the Fourier repres
entation will be very economical:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 41 "ft:=sin(t)+2*cos(5*t)+1/2*sin(sqrt(5)*t);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "plot(ft,t=0..10);" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 16 "fourier(ft,t,w);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 134 "Now, we can't graph the result, but in real life we n
ever want to integrate infinite wave trains. So what happens for a fin
ite signal?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "fw:=int(ft*e
xp(I*w*t),t=-20..20):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "pl
ot([Re(fw),Im(fw)],w=-6..6,color=[red,blue]);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 67 "Clearly the information about the frequency contributi
ons is there." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT 258 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 136 "Extend the \+
time over which the signal is integrated, i.e., consider a longer wave
train, and graph the results. What is your observation?" }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 260 11 "Exercise 3:" }
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 119 "Choose different signals
 and perform the Fourier analysis using the infinite domain, and then \+
different finite domains." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 377 "Note how the ampl
itudes of the signal contributions at different frequencies in the tim
e and frequency domains are preserved. This is a consequence of the li
nearity of the FT. The figure above displays also how we can generate \+
a sequence of oscillatory functions which approximate the Dirac delta \+
(in the limit of an infinite time domain the Dirac delta function is r
ecovered)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 259 26 "Discrete Fourier Trans
form" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 620 "
In the discretized version of the FT we approximate the integral by a \+
Riemann sum. This approximation is not unreasonable, because the integ
rand is assumed to be periodic, which means that the function values a
nd derivatives at the left and right boundaries of the integration int
erval are equal. For this case the Riemann sum is more accurate than s
ome of the usual quadrature rules on an equidistant grid, such as the \+
trapezoid or Simpson rule. The reduction of the infinite domain to a f
inite one means that we are using a Fourier series representation inst
ead of a true Fourier transform [investigate the worksheets " }{TEXT 
19 17 "FourierSeries.mws" }{TEXT -1 2 ", " }{TEXT 19 11 "HeatEqn.mws" 
}{TEXT -1 2 ", " }{TEXT 19 11 "WaveEqn.mws" }{TEXT -1 28 " to find out
 more about it]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 722 "In the discrete FT we not only have a finite domain, but
 in fact we will discretize time and frequency (position and momentum)
. Let us use the position-momentum language familiar from quantum mech
anics, although the same arguments apply in the time-frequency pair of
 variables. If we have a spatial domain of finite extent, then this im
plies that there is a lowest wavenumber (momentum) that can be represe
nted (a sine wave spanning the entire domain with one positive, and on
e negative lobe). Likewise, when we discretize position, there will be
 a maximum wavenumber that can be resolved: a sequence of alternating \+
+1, -1, can be thought of as the highest momentum wave representable f
or a given choice of discretization " }{TEXT 19 2 "dx" }{TEXT -1 103 "
. Let us pick a spacing, and a number of points (which for technical r
easons has to be a power of two)." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "dx:=0.4; n:=8; N:=2^n;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*dx));" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 107 "The spacing in the momentum domain follows from the argu
ments given above: The size of the position domain " }{TEXT 19 4 "N*dx
" }{TEXT -1 167 " controls not only the smallest momentum that can be \+
represented, but also the resolution in the momentum domain. The large
st momentum that can be observed is given by" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 14 "p_max:=N*dp/2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 411 "The momentum and coordinate meshes are arranged in a somewhat \+
unusual way: we begin with the positive domains to be stored in the fi
rst halves of the arrays, followed by the negative values attached in \+
a periodic way (following the largest positive position/momentum come \+
the largest negative positions/momenta). We use vectors from the Linea
rAlgebra package for storage, but lists or arrays would work as well.
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "xv:=Vector(N): pv:=Vect
or(N):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 131 "for i from 1 to \+
N do:\nxv[i]:=evalhf((i-N/2-1)*dx):\nif i < N/2+1 then\npv[i]:=evalhf(
(i-1)*dp): pv[N-i+1]:=evalhf(-i*dp): end if:\nod:" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 7 "pv[20];" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "pointplot(\{seq([i,xv[i]],i=1..N)\});" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "pointplot(\{seq([i,pv[i]],i=1..N)\}
);" }}{PARA 13 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 226 "The meshes include zero as their first point. These are not th
e only possible choices. In particular, one can choose meshes that avo
id the origin, i.e., are spaced a half-step around it, if one has reas
ons to avoid the origin." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 192 "Let us choose a time sequence and generate the tr
ansform. We pick a real function with sine and cosine contributions. T
he function is real-valued, but we prepare for a real and imaginary pa
rt " }{TEXT 19 2 "yR" }{TEXT -1 5 " and " }{TEXT 19 2 "yI" }{TEXT -1 
14 " respectively." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "fx:=s
in(x)+2*cos(2*x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "yR:=Ve
ctor(N): yI:=Vector(N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "for i fr
om 1 to N do: yR[i]:=evalf(subs(x=xv[i],Re(fx))); yI[i]:=evalf(subs(x=
xv[i],Im(fx))); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot
([seq([xv[i],yR[i]],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 37 "It is hard to see the signal! Choose " }{TEXT 19 10 "styl
e=line" }{TEXT -1 17 " to get the idea." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 36 "N_x:=dx*add(yR[i]^2+yI[i]^2,i=1..N);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalhf(FFT(n, var(yR), var(yI)));" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_p:=dp*add(yR[i]^2+yI[i]
^2,i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 68 "This is not normal
ized correctly, so let us define a normalized DFT:" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 124 "yRn:=Vector(N): yIn:=Vector(N): sq:=evalf(
1/sqrt(2*Pi)): for i from 1 to N do: yRn[i]:=dx*sq*yR[i]: yIn[i]:=dx*s
q*yI[i]: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "N_p:=dp*add
(yRn[i]^2+yIn[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 85 "No
w we understand that the FFT function carried out the sum, but did not
 multiply by " }{TEXT 19 2 "dx" }{TEXT -1 79 " (to turn it into the Ri
emann sum), and that it does not include the customary " }{TEXT 19 12 
"1/sqrt(2*Pi)" }{TEXT -1 8 " factor." }}{PARA 0 "" 0 "" {TEXT -1 0 "" 
}}{PARA 0 "" 0 "" {TEXT -1 224 "Observe the efficiency: in one GO the \+
FFT function has calculated N integrals, one for each of the frequency
 (momentum) values represented on the p-mesh. The reason why it is cal
led fast is that the effort does not scale as " }{TEXT 19 3 "N^2" }
{TEXT -1 59 " as would be normally expected from the DFT, but rather a
s " }{TEXT 19 8 "N*log(N)" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 28 "Now let us graph the result:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "plot([seq([pv[i],yRn[i]],i=1
..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "plot
([seq([pv[i],yIn[i]],i=1..N)],style=point,color=blue);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "plot([seq([pv[i],yRn[i]^2+yIn[i]^2]
,i=1..N)],style=point,color=black);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 462 "Evidently, unless the resolution is increased considerably one
 can read off the accurate weight of the different frequency contribut
ions only when a proper lineshape is fitted to the data near the peaks
: the frequency (momentum) grid can easily miss the accurate location \+
of the peak. Nevertheless, we can see in the given example that the si
gnal intensities for the two frequency contributions are approximately
 4:1, as is appropriate for the original function " }{TEXT 19 20 "fx=s
in(x)+2*cos(2*x)" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 344 "Note that the frequency information is r
elatively sharp, because the time interval sampled contains a good num
ber of cycles. The efficiency of the FFT to display information conten
t in terms of population of frequencies is obvious, and this explains \+
the widespread use of FFTs in Electronics, Acoustics, and other branch
es of signal processing." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 450 "There is some redundancy in the above result: We \+
have a degree of symmetry: the real part of the FT is symmetric, while
 the imaginary part is antisymmetric in frequency with the end result \+
of a symmetric intensity profile. The redundancy to some extent is the
 result of using a real signal only. When adding an imaginary part to \+
the function in the time domain (at a different frequency, e.g., f=3),
 the real part of the transform loses its symmetry." }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 261 11 "Exercise 4:" }}{PARA 
0 "" 0 "" {TEXT -1 65 "Play with the above code while keeping the same
 function: change " }{TEXT 19 2 "dx" }{TEXT -1 25 ", while keeping the
 same " }{TEXT 19 1 "N" }{TEXT -1 12 "; then keep " }{TEXT 19 2 "dx" }
{TEXT -1 17 ", while changing " }{TEXT 19 1 "N" }{TEXT -1 203 ". Obser
ve the implications for the information content. If one wishes to work
 with a fixed total length interval, one has to use a window function \+
to cut off the signal outside the desired time interval." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 266 20 "Filtering
 noisy data" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 397 "A classical application of DFT's concerns the removal of noise
 from data. The idea is that noise represents signal components of hig
h frequencies which can be removed by a low-pass filter. In this secti
on we will construct a signal with known smooth components and will ad
d some randomness. Then we will analyze the DFT, remove high frequenci
es and observe the result of an inverse transformation." }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 73 "Maple has a built in
 random number generator for integer randoms. In the " }{TEXT 19 5 "st
ats" }{TEXT -1 82 " package there are also more sophisticated generato
rs, but this one will suffice. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "rand()/10^12;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "evalf(rand()/10^12);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 148 "Subsequent calls result in a sequence of
 pseudorandom numbers. One can changed the seed value for different de
terministic sequences. By dividing by " }{TEXT 19 5 "10^12" }{TEXT -1 
62 " we have uniformly distributed randoms on the interval [0, 1]." }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "fx:=(sin(x/4)+2*cos(1/2*x))
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "yR:=Vector(N): yI:=Vec
tor(N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 118 "for i from 1 to N do: y
R[i]:=evalf(subs(x=xv[i],Re(fx))+1*(rand()/10^12-0.5)); yI[i]:=evalf(s
ubs(x=xv[i],Im(fx))); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
62 "P1:=plot([seq([xv[i],yR[i]],i=1..N)],style=line): display(P1);" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalhf(FFT(n, var(yR), var(
yI)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "plot([seq([pv[i],
yR[i]],i=1..N)],style=line);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 45 "plot([seq([pv[i],yI[i]],i=1..N)],style=line);" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 9 "p_c:=1.5;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 72 "for i from 1 to N do: if abs(pv[i])>p_c then yR[i]:=0
: yI[i]:=0: fi: od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "Now go bac
k to the time domain using the inverse Fourier transform:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "evalhf(iFFT(n, var(yR), var(yI)));
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "P2:=plot([seq([xv[i],yR
[i]],i=1..N)],style=line,color=blue): display(P1,P2);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 122 "Thus, we have demonstrated low-pass filt
ering using a sequence of FFT/iFFT with a truncation step in the frequ
ency domain." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
267 11 "Exercise 5:" }}{PARA 0 "" 0 "" {TEXT -1 83 "Explore the result
s of the above filtering step by changing the cut-off frequency. " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT 262 27 "A quantum mechanics example" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 201 "We proceed with the example of fi
nding the DFT of a Gaussian wavepacket. The example is interesting, be
cause we have a detailed understanding of how the wavepacket is made u
p of plane waves of momenta " }{TEXT 19 1 "p" }{TEXT -1 60 " using a G
aussian weight function defined by a center value " }{TEXT 19 2 "p0" }
{TEXT -1 13 " and a width " }{TEXT 19 2 "dp" }{TEXT -1 1 "." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "x0:=0
; p0:=1/2; sigma:=1/4;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "f
x:=exp(-sigma^2*(x-x0)^2/2)*exp(I*p0*x)/(Pi/(sigma)^2)^(1/4);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "simplify(int(abs(fx)^2,x=-in
finity..infinity));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "The analyt
ical FT to p-space:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "Fp:=
1/sqrt(2*Pi)*int(fx*exp(-I*p*x), x=-infinity..infinity);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 151 "The transform is real-valued, because th
e original function has symmetry in the real part, and anti-symmetry i
n the imaginary part. These combine with " }{TEXT 19 10 "exp(I*p*x)" }
{TEXT -1 154 " in such a way as to cancel the imaginary part of the tr
ansform, since one is integrating an odd function in the imaginary par
t of the combined integrand." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 32 "int(Fp^2,p=-infinity..infinity);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 523 "Now we calculate the DFT. It is advantageous to re-defin
e the position mesh as also cut up in such a way that the array starts
 with small x-values followed by the wrap-around to the large negative
 values. This takes care of the symmetry between the coordinate and mo
mentum representations. If one keeps the mesh used for the time-to-fre
quency domain conversion, an annoying phase (alternating sign) appears
 between adjacent points, and a wrong imaginary part is calculated. Th
is makes it difficult to interpret the results." }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 22 "dx:=0.4; n:=8; N:=2^n;" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*dx));" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 14 "p_max:=N*dp/2;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 29 "xv:=Vector(N): pv:=Vector(N):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 153 "for i from 1 to N do:\nif i < N/2+1 then\nxv[i]:=eva
lhf((i-1)*dx): xv[N-i+1]:=evalhf(-i*dx):\npv[i]:=evalhf((i-1)*dp): pv[
N-i+1]:=evalhf(-i*dp): end if:\nod:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "pv[20];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "
pointplot(\{seq([i,xv[i]],i=1..N)\});" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "pointplot(\{seq([i,pv[i]],i=1..N)\});" }}{PARA 13 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "yR:=Ve
ctor(N): yI:=Vector(N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "for i fr
om 1 to N do: yR[i]:=evalf(subs(x=xv[i],Re(fx))); yI[i]:=evalf(subs(x=
xv[i],Im(fx))); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot
([seq([xv[i],yR[i]],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 57 "plot([seq([xv[i],yI[i]],i=1..N)],style=point,color=
blue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_x:=dx*add(yR[i]
^2+yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "eval
hf(FFT(n, var(yR), var(yI)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 36 "N_p:=dp*add(yR[i]^2+yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 68 "This is not normalized correctly, so let us define a norm
alized DFT:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 124 "yRn:=Vector
(N): yIn:=Vector(N): sq:=evalf(1/sqrt(2*Pi)): for i from 1 to N do: yR
n[i]:=dx*sq*yR[i]: yIn[i]:=dx*sq*yI[i]: od:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 38 "N_p:=dp*add(yRn[i]^2+yIn[i]^2,i=1..N);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "plot([seq([pv[i],yRn[i]],i=1..N)],s
tyle=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "plot([seq([
pv[i],yIn[i]],i=1..N)],style=point,color=blue);" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 118 "The imaginary part is returned as zero, while the r
eal part describes a wavefunction of Gaussian shape centered on p0." }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 263 11 "Exercis
e 6:" }}{PARA 0 "" 0 "" {TEXT -1 103 "Change the parameters of the Gau
ssian wavepacket and observe how they are reflected in the DFT results
." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 215 "As a final check we carry out an error a
nalysis: we do have an analytic answer for the FT, i.e., we can figure
 out how accurate the DFT is for the present example. Let us evaluate \+
the analytic FT at the grid points:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 77 "FT:=Vector(N): for i from 1 to N do: FT[i]:=evalf(sub
s(p= pv[i],Fp), 16); od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "Pick \+
a location on the momentum mesh and compare the results and calculate \+
the relative error:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "i0:=
10: pv[i0];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "FT[i0],yRn[i
0],abs((FT[i0]-yRn[i0])/FT[i0]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
177 "Phenomenal! The agreement is with computational accuracy for the \+
analytic result controlled by the Digits variable; (increase Digits an
d repeat the calculation to confirm this)." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 36 "However, if we go to larger mom
enta:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "i0:=30: pv[i0];" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "FT[i0],yRn[i0],abs((FT[i0]-yRn[i0]
)/FT[i0]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 87 "In fact, if one cal
culates the average relative error over the entire mesh one obtains:" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "add(abs((FT[i]-yRn[i])/FT
[i]),i=1..N)/N;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 221 "A ridiculousl
y large number! This is the result of applying periodic boundary condi
tions in the DFT. Instead of falling to zero the 'Gaussian' as compute
d per DFT connects periodically to the neighbours in adjacent cells." 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "P1:=plot([seq([pv[i],log(
yRn[i])],i=1..N/2)]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "P2
:=plot([seq([pv[i],log(FT[i])],i=1..N/2)],color=blue):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "display([P1,P2],view=[0..4,-25..1])
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 264 11 "Exercise 7:" }}{PARA 0 "" 
0 "" {TEXT -1 112 "Compare the analytic and numerical FT for different
 choices of parameteres, particularly for smaller grid sizes." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 265 11 "Exercise
 8:" }}{PARA 0 "" 0 "" {TEXT -1 117 "Carry out the comparison between \+
analytic and discrete FT for a wavepacket of a different shape, namely
 a Lorentzian." }}{PARA 0 "" 0 "" {TEXT 19 30 "exp(I*p0*x)*a^2/(a^2+(x
-x0)^2)" }{TEXT -1 56 ". Why are the results less accurate than for a \+
Gaussian?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "12
3" 0 }{VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }