{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 }