{VERSION 6 0 "IBM INTEL NT" "6.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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "Blue Emphasis" -1 256 "Times" 0 0 0 0 255 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "Green Emphasis" -1 257 "Times" 1 12 0 128 0 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "Maroon Emphasis" -1 258 "Times" 1 12 128 0 128 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "Dark Red Emphasis" -1 259 "Times" 1 12 128 0 0 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "Red Emphasis" -1 260 "Times" 1 12 255 0 0 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "Purple Emphasis" -1 261 " Times" 1 12 102 0 230 1 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" 260 262 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" 260 263 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" 260 265 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" 260 266 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" 261 267 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 0 0 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 "Heading 1 " -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 128 1 2 1 2 2 2 2 1 1 1 1 } 1 1 0 0 8 4 3 0 3 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times" 1 14 128 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 "Times" 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 }{PSTYLE "" 0 257 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 3 "" 0 "" {TEXT -1 35 "The discrete Fast Fourier Transfo rm" }}{PARA 0 "" 0 "" {TEXT -1 37 "by Peter Stone, Nanaimo, B.C., Cana da" }}{PARA 0 "" 0 "" {TEXT -1 19 "Version: 27.3.2007" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{SECT 1 {PARA 4 " " 0 "" {TEXT -1 60 "Derivation of the discrete Fourier transform and i ts inverse" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{PARA 0 " " 0 "" {TEXT -1 21 "Suppose that we have " }{TEXT 270 1 "n" }{TEXT -1 16 " equally spaced " }{TEXT 268 1 "x" }{TEXT -1 8 " values " } {XPPEDIT 18 0 "x[0],x[1],` . . . `,x[n-1];" "6&&%\"xG6#\"\"!&F$6#\"\" \"%(~.~.~.~G&F$6#,&%\"nGF)F)!\"\"" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "x = n*h" "6#/%\"xG*&%\"nG\"\"\"%\"hGF'" }{TEXT -1 20 ", and corr esponding " }{TEXT 269 1 "y" }{TEXT -1 8 " values " }{XPPEDIT 18 0 "y[ 0],y[1],` . . . `,y[n-1];" "6&&%\"yG6#\"\"!&F$6#\"\"\"%(~.~.~.~G&F$6#, &%\"nGF)F)!\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 81 "We inv estigate the possibilty of interpolating this data by a formula of the form" }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "y(x)=Sum(c[ k]*phi[k](x),k=p..q)" "6#/-%\"yG6#%\"xG-%$SumG6$*&&%\"cG6#%\"kG\"\"\"- &%$phiG6#F/6#F'F0/F/;%\"pG%\"qG" }{TEXT -1 21 " - - - - - - (i), " }}{PARA 0 "" 0 "" {TEXT -1 11 "where each " }{XPPEDIT 18 0 "phi[k](x) " "6#-&%$phiG6#%\"kG6#%\"xG" }{TEXT -1 46 " is a complex exponential f unction of the form" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "phi[k](x)=exp(i*k*alpha*x)/sqrt(n)" "6#/-&%$phiG6#%\"kG6#%\"xG*&-%$ expG6#**%\"iG\"\"\"F(F1%&alphaGF1F*F1F1-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 23 "We choose the constant " }{XPPEDIT 18 0 "alpha" "6#%&alphaG" } {TEXT -1 25 " and the summation range " }{XPPEDIT 18 0 "k = p*` . . `* q" "6#/%\"kG*(%\"pG\"\"\"%&~.~.~GF'%\"qGF'" }{TEXT -1 23 " so that the functions " }{XPPEDIT 18 0 "phi[k](x)" "6#-&%$phiG6#%\"kG6#%\"xG" } {TEXT -1 43 " satisfy a suitable orthogonality relation." }}{PARA 0 " " 0 "" {TEXT -1 4 "Let " }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "S = Sum(phi[k](x[j])*phi[m](x[j]),j = 0 .. n-1);" "6#/%\"SG-%$Su mG6$*&-&%$phiG6#%\"kG6#&%\"xG6#%\"jG\"\"\"-&F+6#%\"mG6#&F06#F2F3/F2;\" \"!,&%\"nGF3F3!\"\"" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 4 "Th en" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 2 " \+ " }{XPPEDIT 18 0 "S = 1/n;" "6#/%\"SG*&\"\"\"F&%\"nG!\"\"" }{TEXT -1 1 " " }{XPPEDIT 18 0 "Sum(exp(k*i*alpha*x[j])*exp(m*i*alpha*x[j]),j = \+ 0 .. n-1);" "6#-%$SumG6$*&-%$expG6#**%\"kG\"\"\"%\"iGF,%&alphaGF,&%\"x G6#%\"jGF,F,-F(6#**%\"mGF,F-F,F.F,&F06#F2F,F,/F2;\"\"!,&%\"nGF,F,!\"\" " }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "`` = 1/n;" "6#/%!G*&\"\"\"F&%\"nG!\"\" " }{TEXT -1 1 " " }{XPPEDIT 18 0 "Sum(exp((k+m)*i*alpha*x[j]),j = 0 .. n-1);" "6#-%$SumG6$-%$expG6#**,&%\"kG\"\"\"%\"mGF,F,%\"iGF,%&alphaGF, &%\"xG6#%\"jGF,/F3;\"\"!,&%\"nGF,F,!\"\"" }{TEXT -1 2 " " }}{PARA 0 " " 0 "" {TEXT -1 2 " " }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "`` = 1/n;" "6#/%!G*&\"\"\"F&%\"nG!\"\"" }{TEXT -1 1 " " } {XPPEDIT 18 0 "Sum(exp((k+m)*i*alpha*j*h),j = 0 .. n-1);" "6#-%$SumG6$ -%$expG6#*,,&%\"kG\"\"\"%\"mGF,F,%\"iGF,%&alphaGF,%\"jGF,%\"hGF,/F0;\" \"!,&%\"nGF,F,!\"\"" }{TEXT -1 3 ", " }}{PARA 0 "" 0 "" {TEXT -1 6 "s ince " }{XPPEDIT 18 0 "x[j]=j*h" "6#/&%\"xG6#%\"jG*&F'\"\"\"%\"hGF)" } {TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 52 "This is a finite geometr ic series with common ratio " }{XPPEDIT 18 0 "r = exp((k+m)*i*alpha*h) ;" "6#/%\"rG-%$expG6#**,&%\"kG\"\"\"%\"mGF+F+%\"iGF+%&alphaGF+%\"hGF+ " }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 5 "Hence" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "S = (1-r^n)/(1-r);" "6#/%\"SG*&,& \"\"\"F')%\"rG%\"nG!\"\"F',&F'F'F)F+F+" }{TEXT -1 2 ", " }}{PARA 0 "" 0 "" {TEXT -1 5 "when " }{XPPEDIT 18 0 "r<>1" "6#0%\"rG\"\"\"" }{TEXT -1 6 ", and " }{XPPEDIT 18 0 "S=1" "6#/%\"SG\"\"\"" }{TEXT -1 6 " when " }{XPPEDIT 18 0 "r=1" "6#/%\"rG\"\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 13 "If we choose " } {XPPEDIT 18 0 "alpha=2*pi/(n*h)" "6#/%&alphaG*(\"\"#\"\"\"%#piGF'*&%\" nGF'%\"hGF'!\"\"" }{TEXT -1 6 ", then" }}{PARA 256 "" 0 "" {TEXT -1 4 " " }{XPPEDIT 18 0 "r = exp((k+m)*2*pi*i/n);" "6#/%\"rG-%$expG6#*,, &%\"kG\"\"\"%\"mGF+F+\"\"#F+%#piGF+%\"iGF+%\"nG!\"\"" }{TEXT -1 2 ". \+ " }}{PARA 0 "" 0 "" {TEXT -1 3 "If " }{XPPEDIT 18 0 "m=-k" "6#/%\"mG,$ %\"kG!\"\"" }{TEXT -1 10 ", we have " }{XPPEDIT 18 0 "r=1" "6#/%\"rG\" \"\"" }{TEXT -1 10 ", so that " }}{PARA 256 "" 0 "" {TEXT -1 3 " " } {XPPEDIT 18 0 "Sum(phi[k](x[j])*phi[m](x[j]),j = 0 .. n-1) = 1;" "6#/- %$SumG6$*&-&%$phiG6#%\"kG6#&%\"xG6#%\"jG\"\"\"-&F*6#%\"mG6#&F/6#F1F2/F 1;\"\"!,&%\"nGF2F2!\"\"F2" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 10 "Otherwise " }{XPPEDIT 18 0 "r^n =1" "6#/)%\"rG%\"nG\"\"\"" } {TEXT -1 5 ", and" }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "Sum(phi[k](x[j])*phi[m](x[j]),j = 0 .. n-1) = 0;" "6#/-%$SumG6$*&-&%$ phiG6#%\"kG6#&%\"xG6#%\"jG\"\"\"-&F*6#%\"mG6#&F/6#F1F2/F1;\"\"!,&%\"nG F2F2!\"\"F<" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 14 "The functions " }{XPPEDIT 18 0 "phi[k](x)=exp(2 *pi*k*x/(n*h))/sqrt(n)" "6#/-&%$phiG6#%\"kG6#%\"xG*&-%$expG6#*,\"\"#\" \"\"%#piGF1F(F1F*F1*&%\"nGF1%\"hGF1!\"\"F1-%%sqrtG6#F4F6" }{TEXT -1 9 " form an " }{TEXT 261 17 "orthonormal basis" }{TEXT -1 33 " for a vec tor space of functions " }{XPPEDIT 18 0 "y(x)" "6#-%\"yG6#%\"xG" } {TEXT -1 52 " of the form (i), with respect to the inner product " }} {PARA 256 "" 0 "" {TEXT -1 3 " f " }{TEXT 264 1 "." }{TEXT -1 2 " g" } {XPPEDIT 18 0 "`` = Sum(f(x[j])*g(x[j]),j = 0 .. n-1);" "6#/%!G-%$SumG 6$*&-%\"fG6#&%\"xG6#%\"jG\"\"\"-%\"gG6#&F-6#F/F0/F/;\"\"!,&%\"nGF0F0! \"\"" }{TEXT -1 4 ", " }}{PARA 0 "" 0 "" {TEXT -1 30 "that is, they \+ are a system of " }{TEXT 261 32 "mutually orthogonal unit vectors" } {TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 45 "These considerations mean that we can expand " }{XPPEDIT 18 0 "y(x)" "6#-%\"yG6#%\"xG" }{TEXT -1 12 " in the form" }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "y(x) = Sum(c[k]*phi[k](x),k = \+ 0 .. n-1);" "6#/-%\"yG6#%\"xG-%$SumG6$*&&%\"cG6#%\"kG\"\"\"-&%$phiG6#F /6#F'F0/F/;\"\"!,&%\"nGF0F0!\"\"" }{TEXT -1 2 ", " }}{PARA 0 "" 0 "" {TEXT -1 6 "where " }{XPPEDIT 18 0 "phi[k](x)=1/sqrt(n)" "6#/-&%$phiG6 #%\"kG6#%\"xG*&\"\"\"F,-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 1 " " } {XPPEDIT 18 0 "exp(2*pi*k*x/(n*h))" "6#-%$expG6#*,\"\"#\"\"\"%#piGF(% \"kGF(%\"xGF(*&%\"nGF(%\"hGF(!\"\"" }{TEXT -1 2 ". " }}{PARA 0 "" 0 " " {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 4 "The " }{TEXT 261 22 "or thogonality relation" }{TEXT -1 1 " " }}{PARA 256 "" 0 "" {TEXT -1 3 " " }{XPPEDIT 18 0 "Sum(phi[k](x[j])*phi[m](x[j]),j = 0 .. n-1) = PIE CEWISE([1, m = -k],[0, m <> -k]);" "6#/-%$SumG6$*&-&%$phiG6#%\"kG6#&% \"xG6#%\"jG\"\"\"-&F*6#%\"mG6#&F/6#F1F2/F1;\"\"!,&%\"nGF2F2!\"\"-%*PIE CEWISEG6$7$F2/F6,$F,F?7$F<0F6,$F,F?" }{TEXT -1 1 " " }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 28 "means that the coeffici ents " }{XPPEDIT 18 0 "c[k]" "6#&%\"cG6#%\"kG" }{TEXT -1 43 " in (i) c an be calculated from the relation" }}{PARA 256 "" 0 "" {TEXT -1 1 " \+ " }{XPPEDIT 18 0 "c[k] = Sum(y[j]*phi[-j](x),j = 0 .. n-1);" "6#/&%\"c G6#%\"kG-%$SumG6$*&&%\"yG6#%\"jG\"\"\"-&%$phiG6#,$F/!\"\"6#%\"xGF0/F/; \"\"!,&%\"nGF0F0F6" }{TEXT -1 1 "," }}{PARA 256 "" 0 "" {TEXT -1 1 " \+ " }{XPPEDIT 18 0 "`` = 1/sqrt(n);" "6#/%!G*&\"\"\"F&-%%sqrtG6#%\"nG!\" \"" }{TEXT -1 1 " " }{XPPEDIT 18 0 "Sum(y[j]*exp(-2*pi*i*j*k/n),j = 0 \+ .. n-1);" "6#-%$SumG6$*&&%\"yG6#%\"jG\"\"\"-%$expG6#,$*.\"\"#F+%#piGF+ %\"iGF+F*F+%\"kGF+%\"nG!\"\"F6F+/F*;\"\"!,&F5F+F+F6" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 5 "Thus " }{XPPEDIT 18 0 "y(x)" "6#-%\"yG6# %\"xG" }{TEXT -1 29 " can be expanded in the form " }}{PARA 256 "" 0 " " {TEXT -1 2 " " }{TEXT 267 2 " " }{XPPEDIT 18 0 "y(x) = 1/sqrt(n); " "6#/-%\"yG6#%\"xG*&\"\"\"F)-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 1 " " } {XPPEDIT 18 0 "Sum(c[k]*exp(2*Pi*i*k*x/(n*h)),k = 0 .. n-1);" "6#-%$Su mG6$*&&%\"cG6#%\"kG\"\"\"-%$expG6#*.\"\"#F+%#PiGF+%\"iGF+F*F+%\"xGF+*& %\"nGF+%\"hGF+!\"\"F+/F*;\"\"!,&F5F+F+F7" }{TEXT -1 16 " ------- (ii) , " }}{PARA 0 "" 0 "" {TEXT -1 6 "where " }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "c[k] = Sum(y[j]*exp(-2*Pi*i*j*k/n),j = 0 .. n -1);" "6#/&%\"cG6#%\"kG-%$SumG6$*&&%\"yG6#%\"jG\"\"\"-%$expG6#,$*.\"\" #F0%#PiGF0%\"iGF0F/F0F'F0%\"nG!\"\"F:F0/F/;\"\"!,&F9F0F0F:" }{TEXT -1 3 ". " }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{TEXT 265 14 "___________ ___" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 17 "Th e coefficients " }{XPPEDIT 18 0 "c[k];" "6#&%\"cG6#%\"kG" }{TEXT -1 17 " constitute the " }{TEXT 261 26 "discrete Fourier transform" } {TEXT -1 15 " of the values " }{XPPEDIT 18 0 "y[0],y[1],` . . . `,y[n- 1];" "6&&%\"yG6#\"\"!&F$6#\"\"\"%(~.~.~.~G&F$6#,&%\"nGF)F)!\"\"" } {TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 7 "Taking " }{XPPEDIT 18 0 "x=x[j]" "6#/%\"xG&F$6#%\"jG" }{XPPEDIT 18 0 "``=j*h" "6#/%!G*&%\"jG\" \"\"%\"hGF'" }{TEXT -1 3 ", " }{XPPEDIT 18 0 "j = 0,` . . . `,n-1;" " 6%/%\"jG\"\"!%(~.~.~.~G,&%\"nG\"\"\"F)!\"\"" }{TEXT -1 29 ", in the re lation (ii) gives " }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "y[j] = Sum(c[k]*exp(2*Pi*i*j*k/n),k = 0 .. n-1);" "6#/&%\"yG6#%\"jG-% $SumG6$*&&%\"cG6#%\"kG\"\"\"-%$expG6#*.\"\"#F0%#PiGF0%\"iGF0F'F0F/F0% \"nG!\"\"F0/F/;\"\"!,&F8F0F0F9" }{TEXT -1 2 ". " }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{TEXT 266 14 "______________" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 38 "which allows the original data values " }{XPPEDIT 18 0 "y[j]" "6#&%\"yG6#%\"jG" }{TEXT -1 39 " \+ to be recovered from the coefficients " }{XPPEDIT 18 0 "c[k];" "6#&%\" cG6#%\"kG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 31 "These equati ons constitute the " }{TEXT 261 34 "inverse discrete Fourier transform " }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 2 " " }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 60 "Formulas for the discrete Fourier transform and its inverse " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 "The " }{TEXT 261 26 "discrete Fourier transform" }{TEXT -1 14 " is defined by" }} {PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "c[k]=1/sqrt(n)" "6#/ &%\"cG6#%\"kG*&\"\"\"F)-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 1 " " } {XPPEDIT 18 0 "Sum(y[j]*exp(-2*Pi*i*j*k/n),j = 0 .. n-1);" "6#-%$SumG6 $*&&%\"yG6#%\"jG\"\"\"-%$expG6#,$*.\"\"#F+%#PiGF+%\"iGF+F*F+%\"kGF+%\" nG!\"\"F6F+/F*;\"\"!,&F5F+F+F6" }{TEXT -1 2 ", " }}{PARA 0 "" 0 "" {TEXT -1 4 "for " }{XPPEDIT 18 0 "j = 1,` . . . `,n" "6%/%\"jG\"\"\"%( ~.~.~.~G%\"nG" }{TEXT -1 5 ", so " }{TEXT 271 1 "n" }{TEXT -1 40 " sam ple points will be transformed into " }{TEXT 272 1 "n" }{TEXT -1 9 " p oints. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 "Let " }{XPPEDIT 18 0 "omega[n] = exp(2*Pi*i/n);" "6#/&%&omegaG6#%\"nG -%$expG6#**\"\"#\"\"\"%#PiGF-%\"iGF-F'!\"\"" }{XPPEDIT 18 0 "`` = cos( 2*pi/n)+i*sin(2*pi/n);" "6#/%!G,&-%$cosG6#*(\"\"#\"\"\"%#piGF+%\"nG!\" \"F+*&%\"iGF+-%$sinG6#*(F*F+F,F+F-F.F+F+" }{TEXT -1 2 ". " }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 "The " }{TEXT 261 18 "nth roots of unity" }{TEXT -1 6 " are: " }{XPPEDIT 18 0 "omega[n]^ 0 = 1,omega[n],omega[n]^2,` . . . `,omega[n]^(n-1);" "6'/*$&%&omegaG6# %\"nG\"\"!\"\"\"&F&6#F(*$&F&6#F(\"\"#%(~.~.~.~G)&F&6#F(,&F(F*F*!\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 19 "They correspond to " } {TEXT 273 1 "n" }{TEXT -1 41 " equally spaced points around the circle " }{XPPEDIT 18 0 "abs(z)=1" "6#/-%$absG6#%\"zG\"\"\"" }{TEXT -1 23 " \+ in the complex plane. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 257 " " 0 "" {TEXT -1 1 " " }{GLPLOT2D 339 330 330 {PLOTDATA 2 "6,-%'CURVESG 6&7S7$$\"\"\"\"\"!$F*F*7$$\"3w\"4hRPij!**!#=$\"3Ikwb#=y_O\"F/7$$\"3E8J #))4-Qn*F/$\"3%[#\\ff*)GLDF/7$$\"3-N5')yke[#*F/$\"3gLj&[K5J!QF/7$$\"3? goz=42`')F/$\"3j9NXR4U7]F/7$$\"3Sb](G._U!zF/$\"3t_H-qceDhF/7$$\"3_\\$R 'oTd#3(F/$\"3!GO#*3qU&fqF/7$$\"3?H$>jubk6'F/$\"3!\\@@&\\!>8\"zF/7$$\"3 VGuIc:x5]F/$\"3-$>p(Qh-a')F/7$$\"37C3nW'R,#QF/$\"3aK+16bcT#*F/7$$\"3$o Y@iF'QDDF/$\"3s&Q\"[M\"oen*F/7$$\"3OB^hAo8X8F/$\"33)fVPO<\"4**F/7$$!3+ qB/u(p5(f!#@$\"2%HhJ<#)******!#<7$$!33)\\T#fB[i8F/$\"3ap%>wGZn!**F/7$$ !3)e:d.:#=YEF/$\"3PA>=\\D`V'*F/7$$!3%*yV1J*3Jx$F/$\"3IzF#e`m3E*F/7$$!3 V(\\AOF:B/&F/$\"3vU'yI2&oN')F/7$$!3[igNw%*[RgF/$\"3!G;)RQ+BqzF/7$$!3/X UwYjJ*3(F/$\"3AfvOg?x_qF/7$$!3&=e_b?/U!zF/$\"3u3HvXQF/7$$!3!p!RS4Nij '*F/$\"3o$p9m#Q%=d#F/7$$!3#p(pT>+.2**F/$\"31V?00]Ug8F/7$$!3/gKG4>&**** *F/$\"3vCA\\O%485$!#?7$$!3__Is=!Q%3**F/$!3q#4*>c=8]8F/7$$!3)>b`sYlSn*F /$!3UNbFGIGKDF/7$$!37_J.8AEj#*F/$!3q4&=j`Bsw$F/7$$!3%\\F)4Kh=u')F/$!3k ENX')3zv\\F/7$$!3EB*z7xng%zF/$!37W(H!pUCrgF/7$$!3XD84Pfc5rF/$!3Y?#o?\" yMJqF/7$$!3!Q%y6Ike[gF/$!3)4j2yeGL'zF/7$$!3g@UD=%)o!*\\F/$!3I_Mh6Mil') F/7$$!3o+1s\"[\"ysPF/$!3#\\IN,%***4E*F/7$$!30yCH\"el'3EF/$!3=V>(fp[Pl* F/7$$!3g67Cvnc\"H\"F/$!3;$RoI*>C;**F/7$$!3EzyOHMQdIFas$!3W1O#>E`*****F /7$$\"3GIAj`Ks(G\"F/$!3lYZ3X=u;**F/7$$\"3EKRqX8/bDF/$!3U$[o%H'z!o'*F/7 $$\"3@%)yb&Q)zNQF/$!36.2<#>x]B*F/7$$\"3O7w4w0I.]F/$!3wHgb5wMe')F/7$$\" 3@\\#)[*R]$4hF/$!3/a*[=H2o\"zF/7$$\"33/wY'Q+$*4(F/$!3)4d)eg?sUqF/7$$\" 3[f=s$Ry,!zF/$!3D1$H)zD(p_v\\F/7$$\"3! G$e@`4+D#*F/$!3i&>2ndo*fQF/7$$\"31mGAp))oa'*F/$!3%)f]nRQ=0EF/7$$\"3@zb 'f`bp!**F/$!3/up91t'4O\"F/7$F($\"36YKhSr8/#)!#F-%'COLOURG6&%$RGBG$\"*+ +++\"!\")$\")AR!)\\F_[lF+-%&STYLEG6#%%LINEG-%'SYMBOLG6#%'CIRCLEG-F$6&7 2F'7$$\"3QnG6D`zQ#*F/$\"3#y*3lBV$o#QF/7$$\"3sva'=\"y1rqF/Fc\\l7$$\"3P) *3lBV$o#QF/F^\\l7$F+F(7$$!3P)*3lBV$o#QF/F^\\l7$$!3sva'=\"y1rqF/Fc\\l7$ $!3QnG6D`zQ#*F/F`\\l7$$!\"\"F*F+7$F`]l$!3#y*3lBV$o#QF/7$F]]lF]]l7$Fj\\ lF`]l7$F+Fc]l7$Ff\\lF`]l7$Fc\\lF]]l7$F^\\lFf]l-%&COLORG6&F\\[l$\"\"%Fd ]lF*$\"\"*Fd]l-Fc[l6#%&POINTG-Fg[l6#%(DIAMONDG-F$6&F\\\\lF^^lFe^l-Fg[l 6#%&CROSSG-F$6&F\\\\lF^^lFe^lFf[l-%%TEXTG6%7$$\"$N\"!\"#$Fd]lFd]lQ*Rea l~axis6\"-F_^l6&F\\[l$F)Fh_lF^`lF^`l-Fc_l6%7$$!#XFh_l$\"$D\"Fh_lQ/Imag inary~axisF[`lF\\`l-%&TITLEG6#%9The~16~th~roots~of~unityG-%*AXESTICKSG 6$\"\"$7$/Fd]l%#-iG/F)%\"iG-%+AXESLABELSG6$Q!F[`lFgal-%%VIEWG6$;$!$D\" Fh_lFf_lF[bl" 1 2 0 1 10 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 0 "Curve 1" "Curve 2" "Curve 3" "Curve 4" "Curve 5" "Curve 6" }}{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 17 "The conjugate of " }{XPPEDIT 18 0 "omega[n] = exp(2*pi/n)" "6#/&%&omegaG6#%\"nG-%$expG6#*(\"\"#\"\" \"%#piGF-F'!\"\"" }{TEXT -1 4 " is " }{XPPEDIT 18 0 "conjugate(omega[n ]) = exp(-2*pi/n);" "6#/-%*conjugateG6#&%&omegaG6#%\"nG-%$expG6#,$*(\" \"#\"\"\"%#piGF1F*!\"\"F3" }{TEXT -1 13 ", and we have" }}{PARA 256 " " 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "c[k]=1/sqrt(n)" "6#/&%\"cG6#%\" kG*&\"\"\"F)-%%sqrtG6#%\"nG!\"\"" }{XPPEDIT 18 0 "Sum(y[j]*conjugate(o mega[n])^(j*k),j = 0 .. n-1);" "6#-%$SumG6$*&&%\"yG6#%\"jG\"\"\")-%*co njugateG6#&%&omegaG6#%\"nG*&F*F+%\"kGF+F+/F*;\"\"!,&F3F+F+!\"\"" } {TEXT -1 2 ". " }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{TEXT 262 14 "____ __________" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 42 "In a similar way, the inverse transform is" }} {PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "y[j] = 1/sqrt(n);" "6 #/&%\"yG6#%\"jG*&\"\"\"F)-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 2 " " } {XPPEDIT 18 0 "Sum(c[k]*omega[n]^(j*k),k = 0 .. n-1);" "6#-%$SumG6$*&& %\"cG6#%\"kG\"\"\")&%&omegaG6#%\"nG*&%\"jGF+F*F+F+/F*;\"\"!,&F0F+F+!\" \"" }{TEXT -1 2 ". " }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{TEXT 263 14 "______________" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT 261 17 "Other conventions" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 14 "Interchanging " }{XPPEDIT 18 0 "omega[n]" "6#&% &omegaG6#%\"nG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "conjugate(omega[n] );" "6#-%*conjugateG6#&%&omegaG6#%\"nG" }{TEXT -1 97 " in these two fo rmulas just reverses the (cyclical) order of the values in each of the sequences " }{XPPEDIT 18 0 "c[k]" "6#&%\"cG6#%\"kG" }{TEXT -1 9 " and the " }{XPPEDIT 18 0 "y[j]" "6#&%\"yG6#%\"jG" }{TEXT -1 172 ". This l eads to the convention commonly used in the physical sciences, as oppo sed to the original choice, which is the convention generally used in \+ electrical engineering. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 47 "Another common variation is to omit the factor " } {XPPEDIT 18 0 "1/sqrt(n)" "6#*&\"\"\"F$-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 35 " in the first formula and replace " }{XPPEDIT 18 0 "1/sqrt(n) " "6#*&\"\"\"F$-%%sqrtG6#%\"nG!\"\"" }{TEXT -1 4 " by " }{XPPEDIT 18 0 "1/n" "6#*&\"\"\"F$%\"nG!\"\"" }{TEXT -1 43 " in the formula for th e inverse transform." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 36 "The fast Fourier transform algorithm" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{PARA 0 "" 0 "" {TEXT -1 43 "Suppose that we have n complex data values " }{XPPEDIT 18 0 "a[ 0],a[1],` . . . `,a[n-1];" "6&&%\"aG6#\"\"!&F$6#\"\"\"%(~.~.~.~G&F$6#, &%\"nGF)F)!\"\"" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 23 "The t ransformed points " }{XPPEDIT 18 0 "A[0],A[1],` . . . `,A[n-1];" "6&&% \"AG6#\"\"!&F$6#\"\"\"%(~.~.~.~G&F$6#,&%\"nGF)F)!\"\"" }{TEXT -1 17 ", are defined by " }}{PARA 256 "" 0 "" {TEXT -1 5 " " }{XPPEDIT 18 0 "A[j] = Sum(a[k]*exp(2*Pi*i*j*k/n),k = 0 .. n-1);" "6#/&%\"AG6#%\"jG -%$SumG6$*&&%\"aG6#%\"kG\"\"\"-%$expG6#*.\"\"#F0%#PiGF0%\"iGF0F'F0F/F0 %\"nG!\"\"F0/F/;\"\"!,&F8F0F0F9" }{TEXT -1 2 ", " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "`` = \+ Sum(a[k]*omega[n]^(j*k),k = 0 .. n-1);" "6#/%!G-%$SumG6$*&&%\"aG6#%\"k G\"\"\")&%&omegaG6#%\"nG*&%\"jGF-F,F-F-/F,;\"\"!,&F2F-F-!\"\"" }{TEXT -1 13 " ------- (i)," }}{PARA 0 "" 0 "" {TEXT -1 6 "where " }{XPPEDIT 18 0 "omega[n] = exp(2*Pi*i)/n" "6#/&%&omegaG6#%\"nG*&-%$expG6#*(\"\"# \"\"\"%#PiGF.%\"iGF.F.F'!\"\"" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 43 "Direct computation of all the coefficients " }{XPPEDIT 18 0 "A[j]" "6#&%\"AG6#%\"jG" }{TEXT -1 29 " requires the calculation \+ of " }{XPPEDIT 18 0 "n^2" "6#*$%\"nG\"\"#" }{TEXT -1 19 " terms of the form " }{XPPEDIT 18 0 "a[k]*omega[n]^(j*k)" "6#*&&%\"aG6#%\"kG\"\"\") &%&omegaG6#%\"nG*&%\"jGF(F'F(F(" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 14 "Assuming that " }{TEXT 275 1 "n" }{TEXT -1 62 " is even, \+ the sum (i) can be split into two sums based on the " }{TEXT 260 4 "ev en" }{TEXT -1 5 " and " }{TEXT 260 3 "odd" }{TEXT -1 21 " values of th e index " }{TEXT 274 1 "k" }{TEXT -1 1 "." }}{PARA 256 "" 0 "" {TEXT -1 2 " " }{XPPEDIT 18 0 "A[j]=Sum(a[k]*omega[n]^(j*k),k = 0 .. n-1)" "6#/&%\"AG6#%\"jG-%$SumG6$*&&%\"aG6#%\"kG\"\"\")&%&omegaG6#%\"nG*&F'F0 F/F0F0/F/;\"\"!,&F5F0F0!\"\"" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "`` = Sum(a[ 2*r]*omega[n]^(j*``(2*r)),r = 0 .. n/2-1)+Sum(a[2*r+1]*omega[n]^(j*``( 2*r+1)),r = 0 .. n/2-1);" "6#/%!G,&-%$SumG6$*&&%\"aG6#*&\"\"#\"\"\"%\" rGF/F/)&%&omegaG6#%\"nG*&%\"jGF/-F$6#*&F.F/F0F/F/F//F0;\"\"!,&*&F5F/F. !\"\"F/F/F@F/-F'6$*&&F+6#,&*&F.F/F0F/F/F/F/F/)&F36#F5*&F7F/-F$6#,&*&F. F/F0F/F/F/F/F/F//F0;F=,&*&F5F/F.F@F/F/F@F/" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "`` = Sum(a[2*r]*omega[n]^(j*``(2*r)),r = 0 .. n/2-1)+omega[n]^j*Sum (a[2*r+1]*omega[n]^(j*``(2*r)),r = 0 .. n/2-1);" "6#/%!G,&-%$SumG6$*&& %\"aG6#*&\"\"#\"\"\"%\"rGF/F/)&%&omegaG6#%\"nG*&%\"jGF/-F$6#*&F.F/F0F/ F/F//F0;\"\"!,&*&F5F/F.!\"\"F/F/F@F/*&)&F36#F5F7F/-F'6$*&&F+6#,&*&F.F/ F0F/F/F/F/F/)&F36#F5*&F7F/-F$6#*&F.F/F0F/F/F//F0;F=,&*&F5F/F.F@F/F/F@F /F/" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "`` = A[j]^e+omega[n]^j*A[j]^o;" "6#/ %!G,&)&%\"AG6#%\"jG%\"eG\"\"\"*&)&%&omegaG6#%\"nGF*F,)&F(6#F*%\"oGF,F, " }{TEXT -1 14 " ------- (ii)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "A[j]^e" "6#)&%\"AG6#%\" jG%\"eG" }{TEXT -1 13 " denotes the " }{TEXT 276 1 "j" }{TEXT -1 49 " \+ th component of the Fourier transform of length " }{XPPEDIT 18 0 "n/2 " "6#*&%\"nG\"\"\"\"\"#!\"\"" }{TEXT -1 18 ", formed from the " } {TEXT 260 4 "even" }{TEXT -1 28 " components of the original " } {XPPEDIT 18 0 "a[k]" "6#&%\"aG6#%\"kG" }{TEXT -1 4 "'s. " }}{PARA 0 " " 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "A[j]^o" "6#)&%\"AG6#%\"jG%\"oG" }{TEXT -1 13 " denotes the " }{TEXT 277 1 "j" }{TEXT -1 49 " th compon ent of the Fourier transform of length " }{XPPEDIT 18 0 "n/2" "6#*&%\" nG\"\"\"\"\"#!\"\"" }{TEXT -1 18 ", formed from the " }{TEXT 260 3 "od d" }{TEXT -1 28 " components of the original " }{XPPEDIT 18 0 "a[k]" " 6#&%\"aG6#%\"kG" }{TEXT -1 3 "'s." }}{PARA 0 "" 0 "" {TEXT -1 58 "Each of these transforms would require the calculation of " }{XPPEDIT 18 0 "(n/2)^2=n^2/4" "6#/*$*&%\"nG\"\"\"\"\"#!\"\"F(*&F&F(\"\"%F)" } {TEXT -1 19 " terms of the form " }{XPPEDIT 18 0 "a[k]*omega[n]^(j*k) " "6#*&&%\"aG6#%\"kG\"\"\")&%&omegaG6#%\"nG*&%\"jGF(F'F(F(" }{TEXT -1 17 ", for a total of " }{XPPEDIT 18 0 "n^2/2" "6#*&%\"nG\"\"#F%!\"\"" }{TEXT -1 42 " such computations. Thus calculating the " }{XPPEDIT 18 0 "A[j]" "6#&%\"AG6#%\"jG" }{TEXT -1 47 " from these two transforms using (ii) requires " }{TEXT 261 31 "half the number of calculations " }{TEXT -1 29 " that are required using (i)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 115 "This process can be repe ated eventually arriving at an algorithm for which the calculation tim e is proportional to " }{XPPEDIT 18 0 "n*log[2](n)" "6#*&%\"nG\"\"\"-& %$logG6#\"\"#6#F$F%" }{TEXT -1 13 " rather than " }{XPPEDIT 18 0 "n^2 " "6#*$%\"nG\"\"#" }{TEXT -1 108 ". The fast Fourier transform (FFT) a lgorithm can make a huge saving in computation time for large values o f " }{TEXT 278 1 "n" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 44 "First versions of the fast Four ier transform" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "This version is similar to a version which is in the Mapl e library." }}{PARA 0 "" 0 "" {TEXT -1 62 "It uses lists for input and output but uses arrays internally." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 822 "fft := proc(B,m)\nlocal A, n,l,le,le1,u,w,j,i,ip,t,s,nv2,nm1,k,pl;\n A := convert(B,array);\n \+ n := 2^m; le1 := n;\n for l from 1 to m do\n le := le1;\n \+ le1 := 1/2*le;\n u := 1.0; \n pl := Pi/le1;\n w := eva lf(cos(pl)) - I*evalf(sin(pl));\n for j from 1 to le1 do\n \+ for i from j to n by le do\n ip := i+le1;\n t \+ := A[i]+A[ip];\n s := A[i]-A[ip]; \n A[ip] := s* u;\n A[i] := t;\n end do;\n u := u*w;\n \+ end do;\n end do;\n nv2 := 1/2*n;\n nm1 := n-1;\n j := 1;\n for i from 1 to nm1 do\n if i " 0 "" {MPLTEXT 1 0 821 "ifft := proc(B,m)\nlocal A,n,l,le,le1,u,w,j,i,ip,t,s,nv2,nm1,k;\n A := convert(B,array);\n n := 2^m; le1 := n;\n for l from 1 \+ to m do\n le := le1; le1 := 1/2*le; u := 1.0;\n w := eva lf(cos(Pi/le1)) + I*evalf(sin(Pi/le1));\n for j from 1 to le1 d o\n for i from j to n by le do\n ip := i+le1 ;\n t := A[i]+A[ip]; s := A[i]-A[ip]; \n \+ A[ip] := s*u; A[i] := t;\n end do;\n u := u*w; \n end do;\n end do;\n nv2 := 1/2*n; nm1 := n-1; j := 1; \n for i from 1 to nm1 do\n if i " 0 "" {MPLTEXT 1 0 74 "interface(verboseproc=2):\neval(FFT );\neval(iFFT);\ninterface(verboseproc=1):" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 8 "Examples" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 9 "Example 1" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 414 "# set up the data\nn := 4;\ny := array(0..n-1,[7.,5.,6.,9.]):\n convert(y,list);\n\n# construct transformed data\nC := array(0..n-1): \nw := evalf(exp(-2*Pi*I/n));\nj := 'j':\nfor k from 0 to n-1 do\n C [k] := sum(y[j]*w^(j*k),j=0..n-1);\nend do:\nevalf(convert(C,list));\n \n# construct inverse transform\ny := array(0..n-1):\nk := 'k':\nfor j from 0 to n-1 do\n y[j]:=sum(C[k]*w^(-j*k),k=0..n-1)/n;\nend do:\ne valf(convert(y,list));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\"% " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7&$\"\"(\"\"!$\"\"&F&$\"\"'F&$\"\" *F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"wG^#$!\"\"\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7&$\"#F\"\"!^$$\"\"\"F&$\"\"%F&$!\"\"F&^$F($ !\"%F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7&^$$\"+++++q!\"*$\"\"!F)^$$ \"+++++]F'F(^$$\"+++++gF'F(^$$\"+++++!*F'F(" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "fft([7.,5.,6.,9.], 2);\nifft(%,2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7&$\"#F\"\"!^$$\"#5 !\"\"$\"#SF*$!#5F*^$$\"$+\"!\"#$!$+%F2" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7&^$$\"+++++q!\"*$\"\"!F)^$$\"+++++]F'F(^$$\"+++++gF'$!\"!F)^$$ \"+++++!*F'F0" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 9 "Exampl e 2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 358 "# set up the data\nn := 8;\nf := x -> sin(x+0.3)+I*c os(2*x+1.)/2+sin(x/2+0.2)/3+1;\nh := evalf(Pi/4);\nxvals :=[seq(i*h,i= 0..n-1)]:\nyvals := map(f,xvals);\ny := array(0..n-1,yvals):\n\n# cons truct transformed data\nC := array(0..n-1):\nw := evalf(exp(-2*Pi*I/n) );\nj := 'j':\nfor k from 0 to n-1 do\n C[k] := sum(y[j]*w^(j*k),j=0 ..n-1);\nend do:\nevalf(convert(C,list));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\")" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGf* 6#%\"xG6\"6$%)operatorG%&arrowGF(,*-%$sinG6#,&9$\"\"\"$\"\"$!\"\"F2F2* &^##F2\"\"#F2-%$cosG6#,&F1F9$F2\"\"!F2F2F2*&#F2F4F2-F.6#,&F1F8$F9F5F2F 2F2F2F2F(F(F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+N;)R&y!#5" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%&yvalsG7*^$$\"+RF,$!+I\\N2UF,^$$\"+&31%)G #F,$!+K:^,FF,^$$\"+.nvsfF,$\"+I\\N2UF," }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"wG^$$\"+5y1rq!#5$!+5y1rqF(" }}{PARA 12 "" 1 "" {XPPMATH 20 " 6#7*^$$\"+9*)f3(*!\"*$!\"&!#5^$$\"+;Dy:nF*$!+FROPSF'^$$!+8y._ " 0 "" {MPLTEXT 1 0 32 "yv als;\nfft(yvals,3);\nifft(%,3);\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#7 *^$$\"+ RF*$!+I\\N2UF*^$$\"+&31%)G#F*$!+K:^,FF*^$$\"+.nvsfF*$\"+I\\N2UF*" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#7*^$$\"+9*)f3(*!\"*$!\"&!#5^$$\"+0Dy:n F*$!+FROPSF'^$$!+=y._*G/\"F?$\"*It5(GF*^$$\"+Ch%Qh\"F'$\"+*)RG_6F'^$$\"+2Dy: nF*$\"+GROPSF'" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#7*^$$\"+=Luh8!\"*$\" +H:^,F!#5^$$\"+T)*oq?F'$!+E\\N2UF*^$$\"+1s;LAF'$!+H:^,FF*^$$\"+0Us$z\" F'$\"+C\\N2UF*^$$\"+_'o6.\"F'F(^$$\"+;T*)>RF*$!+C\\N2UF*^$$\"+&41%)G#F *$!+M:^,FF*^$$\"+4nvsfF*$\"+E\\N2UF*" }}}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{SECT 1 {PARA 4 " " 0 "" {TEXT -1 9 "Example 3" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 320 "# set up the data\nn := 8; \nf := x -> 1/(1+x^2+I*x);\nh := 0.1;\nxvals :=[seq(i*h,i=0..n-1)]:\ny vals := map(f,xvals);\ny := array(0..n-1,yvals):\n\n# construct transf ormed data\nC := array(0..n-1):\nw := evalf(exp(-2*Pi*I/n));\nj := 'j' :\nfor k from 0 to n-1 do\n C[k] := sum(y[j]*w^(j*k),j=0..n-1);\nend do:\nevalf(convert(C,list));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"n G\"\")" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGf*6#%\"xG6\"6$%)operat orG%&arrowGF(*&\"\"\"F-,(F-F-*$)9$\"\"#F-F-*&F1F-^#F-F-F-!\"\"F(F(F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"\"\"!\"\"" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%&yvalsG7*^$$\"+++++5!\"*$!\"!\"\"!^$$\"+8L([!)*!#5 $!+g`z2(*!#6^$$\"+.zYs#*F0$!+/p;$y\"F0^$$\"+sTGG&)F0$!+wTBZBF0^$$\"+2' pXq(F0$!+S\"[nl#F0^$$\"+C0$He#F0" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"wG^$$ \"+5y1rq!#5$!+5y1rqF(" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#7*^$$\"+)*e'f Q'!\"*$!+e.\\\"e\"F'^$$\"+Z(pAd%!#5$!+\\(*[iNF-^$$\"+1g!zZ$F-$!*1kML)F -^$$\"+I_I!)GF-$\"*j9Yx%F-^$$\"+'ROVS#F-$\"+giA/:F-^$$\"+![/J)=F-$\"+g 0(4k#F-^$$\"+7rPw5F-$\"+c6. " 0 "" {MPLTEXT 1 0 32 "yvals; \nfft(yvals,3);\nifft(%,3);\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#7*^$$ \"+++++5!\"*$!\"!\"\"!^$$\"+8L([!)*!#5$!+g`z2(*!#6^$$\"+.zYs#*F.$!+/p; $y\"F.^$$\"+sTGG&)F.$!+wTBZBF.^$$\"+2'pXq(F.$!+S\"[nl#F.^$$\"+C0$He#F." }} {PARA 12 "" 1 "" {XPPMATH 20 "6#7*^$$\"+,f'fQ'!\"*$!+d.\\\"e\"F'^$$\"+ W(pAd%!#5$!+N(*[iNF-^$$\"+**f!zZ$F-$!*!RYL$)F-^$$\"+7_I!)GF-$\"*x9Yx%F -^$$\"+qjL/CF-$\"+riA/:F-^$$\"+WW5$)=F-$\"+h0(4k#F-^$$\"+hqPw5F-$\"+I6 .^$$\"+9L([!)*!#5$!+p`z2(*!#6^$$\"+/zYs#*F.$!+ 0p;$y\"F.^$$\"+rTGG&)F.$!+wTBZBF.^$$\"+4'pXq(F.$!+T\"[nl#F.^$$\"+E0$He#F." }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}{SECT 1 {PARA 0 "" 0 "" {TEXT -1 24 "Cod e for drawing picture" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" } }}{SECT 1 {PARA 0 "" 0 "" {TEXT -1 20 "16 th roots of unity" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 442 "p ts := [seq([cos(i*Pi/8),sin(i*Pi/8)],i=0..15)]:\np1 := plot([[cos(t),s in(t),t=0..2*Pi],pts$3],style=[line,point$3],\n symbol=[circle,di amond,cross],color=[coral,COLOR(RGB,.4,0,.9)$3]):\nt1 := plots[textplo t]([[1.35,-0.1,`Real axis`],\n [-.45,1.25,`Imaginary axis`]], color=COLOR(RGB,.01,.01,.01)):\nplots[display]([p1,t1],view=[-1.25..1. 35,-1.25..1.35],\n title=`The 16 th roots of unity`,tickmarks=[3,3],y tickmarks=[-1=`-i`,1=`i`]);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 ";" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "4 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }