习题16.4Fourier变换和Fourier积分1.求下列定义在(-o,+o)的函数的Fourier变换：A 00;{o，其它;Te-2x,x≥0,(3) f(x)=e-ax, α>0;(4)f(x)=0,x0o[Aeo*dx=(1-e-o)解(1) (o)=[ f(x)e-*dx=]io(2) F(o)= ft f(x)e-odx = ft e-(ao)dx+ [e(a-i0)dx20a+io a-io"a+?(3) (o)= f(x)e-iox dx = [ e-ar-ioxdx= Je-arcoswxdxot.dt（利用例15.2.8的结果）cOsVa"Va)ENaVa(4) F(0)= J f(x)e-ior dx= J" e-(2+o)dx = 2+i0(5) F(0)=f f(x)e*dx=f, Acos 0xe-o* dxAcosxcoswxdx（虚部为奇函数，积分为0）,[cos(o -0)x + cos(0 +0)x]dxsin(0-0.)8, sin(0+0)8(0-0)(0+0)2.求f(x)=e-a（xe[0,+)，α>0）的正弦变换和余弦变换。解正弦变换：F(o)= J, (x)sin oxdx = f, e" sin oxdx =α2+0?余弦变换：
习题 16.4 Fourier 变换和 Fourier 积分 1．求下列定义在(−∞,+∞) 的函数的 Fourier 变换： ⑴ ⎩ ⎨ ⎧ 0 ; ⑶ f x a x ( ) = e− 2 ， a > 0 ; ⑷ ⎩ ⎨ ⎧ ≤ = 0, | | ; cos , | | , ( ) 0 δ ω δ x A x x f x ω0 ≠ 0是常数， ω0 π δ = 。 解 （1） ( ) ( ) i x f f x e dx ω ω +∞ − −∞ = ∫ 0 i x Ae dx δ − ω = ∫ = (1 ) ωδ ω i e i A − − 。 （2） ( ) ( ) i x f f x e dx ω ω +∞ − −∞ = ∫ 0 ( ) ( ) 0 a i x a i x e dx e ω ω +∞ − + − −∞ = + ∫ ∫ dx 1 1 a iω a iω = + + − = 2 2 2 a +ω a 。 （3） ( ) ( ) i x f f x e dx ω ω +∞ − −∞ = ∫ 2 ax i x e dx ω +∞ − − −∞ = = ∫ 2 cos ax e x ω +∞ − ∫−∞ dx 2 0 2 cos t t t e d a a +∞ − ω = ∫ (利用例 15.2.8 的结果) 2 2 a e a ω π ⎛ ⎞ −⎜ ⎟ ⎝ ⎠ = = a e a 4 2 ω π − 。 （4） ( ) ( ) i x f f x e dx ω ω +∞ − −∞ = ∫ (2 ) 0 i x e dx ω +∞ − + = = ∫ 2 + iω 1 。 （5） ( ) ( ) i x f f x e dx ω ω +∞ − −∞ = ∫ = 0 cos i x A xe dx δ ω δ ω − ∫− 0 Acos x cos x δ δ ω ω − = ∫ dx（虚部为奇函数，积分为 0） 0 0 [cos( ) cos( ) ] 2 A x x dx δ δ ω ω ω ω − = − + + ∫ ＝ 0 0 0 0 sin( ) sin( ) ( ) ( ) A ω ω δ ω ω δ ω ω ω ω ⎡ − + ⎤ ⎢ + ⎥ ⎣ − + ⎦ 。 2．求 f x( ) = e− ax（ x ∈[0,+∞)，a > 0）的正弦变换和余弦变换。 解 正弦变换： 0 f ( ) ω ω f x( )sin xdx +∞ = ∫ 0 sin ax e xdx 2 2 ω ω a + ω , +∞ − = = ∫ 余弦变换： 1
aF(o)= f f(x)cos wxdx = f e" cos oxdx =α?+0?元Je, x≥0,Jsinx,0≤x≤3. 设()=[4时，(x-t)=e-(-)，所以2Tre sin()d =le(+)F(x)=e-x2[e-(x-), x>1,当0le-*(1+e?),222
0 f ( ) ω ω f x( ) cos xdx +∞ = ∫ 0 cos ax e xdx 2 2 a +ω a ω 。 +∞ − = = ∫ 3．设 ⎩ ⎨ ⎧ 时， f x 1( ) − = t e− − ( x t) ，所以 2 2 0 1 ( ) sin( ) (1 ) 2 x t x F x e e t dt e e π π − − = = ∫ + ； 当0 2 x π − = ⎨ ⎩ ≤ ，所以 0 1 ( ) sin( ) (sin cos ) 2 x x t x F x e e t dt x x e − − = = − ∫ + 。 于是 ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ + > − + < ≤ ≤ ∗ = − − . 2 (1 ), 2 1 , 2 (sin cos ), 0 2 1 0, 0, ( ) 2 1 2 π π π e e x x x e x x f f x x x 2
习题16.5快速Fourier变换N-I1．说明离散Fourier变换X(i)=2x(n)e/%可以看成Fourier 变换f(o)-f f(x)e-o dx的离散近似形式的推广。解假设の>0-i(o)=[t f(x)e_ dx ~Sf(nAr)eAr7==002元i取A使，则k为整数时，w+n=W"。于是记W=eNN2元(o)=w"f(kN +n)Ar)Ax,n=0k=of(kN +n)Ar)Ax，所以记x(n)=ka-Ne-2m%Zwx(n) =)X()= f(jo) =n=0n=02.证明正交关系式2/hk1-2元%Ze=oikNO-2.nk2元/nkI解显然，j=k时，Ne=1YZ0下面考虑j+k，不妨设j<k。根据当+1是方程x=1的一个根时，2元;(k-)N-1令=e则=N=e2(k-)x=1。于是有"=0，去1，#m=0*-NNe=09e"ZNLn=03.设N=pq（p,qeN），构造只需o(p+q)N)次运算的Fourier变换3
习题 16.5 快速 Fourier 变换 1． 说明离散 Fourier 变换 X j x n i n j N n N ( ) = ( ) e − = − ∑ 2 0 1 π 可以看成 Fourier 变换 ∫ +∞ −∞ − f = f x dx iωx (ω) ( )e ˆ 的离散近似形式的推广。 解 假设ω > 0 2 ˆ 2 ( ) ( ) e x i f f x dx ω π ω π +∞ − −∞ = ∫ 2 ( ) 2 ( ) e n x i n f n x x ω π π +∞ ∆ − =−∞ ≈ ∑ ∆ ∆ ， 取∆x 使 1 2 x N ω π ∆ = ，记 2 e i W N π − = ，则k 为整数时，W kN +n = W n 。于是 1 0 ˆ( ) (( ) ) N n n k f ω W f kN n x − +∞ = =−∞ = + ∑ ∑ ∆ ∆x ， 记 ( ) (( ) ) k x n f kN n x +∞ =−∞ = ∑ + ∆ ∆x，所以 ˆ X ( )j f = ( jω) 1 0 ( ) N jn n W x n − = = ∑ 1 2 0 ( ) e N nj i N n x n π − − = = ∑ 。 2． 证明正交关系式 j k N nk i N n N n j i N , 2 1 0 2 e e 1 δ π π ∑ = − = − 。 解 显然， j = k 时， 1 2 2 0 1 e e N nk nk i i N N N n π π − − = ∑ =1。 下面考虑 j ≠ k ，不妨设 j < k 。根据当ξ ≠ 1是方程 的一个根时， 有 ，令 = 1 N x 0 1 0 ∑ = − = N n n ξ ( ) 2 e 1 k j i N π ξ − = ≠ ，则 2( ) e N k j π i ξ − = =1。于是 1 0 N n n ξ − = ∑ = 1 1 ( ) 2 2 2 0 0 1 1 e e e N N n j nk n k j i i i N N N N N n n π π π − − − − = = ∑ ∑ = = 0。 3． 设 N = pq（ p, q ∈N），构造只需O(( p + q)N) 次运算的 Fourier 变换 3
算法。2元1解令W=e，则k为整数时，wk+"=W"。假设j= jq+Jo, Jj, =0,1,.,p-1, J=0,1,.,q-1,n=np+no,n, =0,1,..-,q-1, no=0,1,",p-1。N!Zx(n)wnEx(n)e-2%X(i)=Nn=0h=0nZx(np+ng)wj(mp+n)20=0m=0FwmgEx(np+n)wnip。m=0no=0固定」，计算艺x(mp+n)wmp需要q-1次乘法（n=0不需要做乘n=0x(mpP+n)Wp是相同的，无需重复计算，所法），对于相同的j。，420有此类和式共需q(q-1)次乘法。对n.求和需要(p-1)N次乘法，所以，总共需要q(q-1)+(p-1)N=O(p+q)N)次乘法。4．对N=23，具体写出以2为底的FFT的计算流程。2元ini解记W=e=e-，则w4=-1,W=l。可得计算公式X()-Zx(n)W", j=,1,..,7.n=0=[x(0) +(-1) x(4)) + W' [x(1) +(-1) x(5)]+W2j ([x(2)+(-1) x(6))+ W'[x(3)+ (-1) x(7)]) 。计算流程第一步：x(i) = x(i)+ x(i+ 4),x(i+ 4)=W'[x(i)-x(i+ 4)],i=0,1,2,34
算法。 解 令 2 e i W N π − = ，则k 为整数时，W kN +n = W n 。假设 1 0 1 0 j j = +q j , 0 j = ,1," " , p −1, j = 0,1, , q −1, 1 0 1 0 n n = + p n , n = 0,1," " , q −1, n = 0,1, , p −1。 X j x n i n j N n N ( ) = ( ) e − = − ∑ 2 0 1 π 1 0 ( ) N jn n x n W − = = ∑ 1 0 0 1 1 1 ( ) 1 0 0 0 ( ) p q j n p n n n x n p n W − − + = = = + ∑ ∑ 0 1 0 1 1 1 1 0 0 0 ( ) p q jn n j p n n W x n p n W − − = = = + ∑ ∑ 0 。 固定 j ，计算 1 0 需要 1 1 1 0 0 ( ) q n j p n x n p n W − = ∑ + q −1次乘法（ 1 n = 0不需要做乘 法），对于相同的 ， 是相同的，无需重复计算，所 有此类和式共需 次乘法。对 求和需要 0 j 1 0 1 1 1 0 0 ( ) q n j p n x n p n W − = ∑ + q q( −1) 0 n ( p −1)N 次乘法，所以， 总共需要q q( 1− +) ( p −1)N = O(( p + q)N) 次乘法。 4． 对 N = 23，具体写出以 2 为底的 FFT 的计算流程。 解 记 2 8 4 e e i i W π π − − = = ，则 4 8 W = −1,W =1。可得计算公式 7 0 ( ) ( ) , 0,1, ,7. jn n X j x n W j = = = ∑ " [ (0) ( 1) (4)] [ (1) ( 1) (5)] j j j = + x x − +W x + − x 2 {[ (2) ( 1) (6)] [ (3) ( 1) (7)]} j j j j + + W x − x +W x + − x 。 . 计算流程 第一步： 1 1 ( ) ( ) ( 4), ( 4) [ ( ) ( 4)], 0,1, 2,3 i x i x i x i x i W x i x i i = + + + = − + = 4
第二步：x(i) = x(i)+ x,(i+ 2),x(i+2)= W2[x,(i)-x(i+2)],i= 0,1,4,5第三步：X(i)= x2(i)+ x2(i+1),X(i+4)= x(i)-x(i+1),i=0,2,X(i)= x (i+3)+ x(i+4),X(i+ 4) = x (i+3) -x,(i+ 4),i = 1,35
第二步： 2 1 1 2 2 1 1 ( ) ( ) ( 2), ( 2) [ ( ) ( 2)], 0,1, 4,5 i x i x i x i x i W x i x i i . = + + + = − + = 第三步： 2 2 2 2 2 2 2 2 ( ) ( ) ( 1), ( 4) ( ) ( 1), 0, 2, ( ) ( 3) ( 4), ( 4) ( 3) ( 4), 1,3 X i x i x i X i x i x i i X i x i x i X i x i x i i . = + + + = − + = = + + + + = + − + = 5
计算实习题（在教师的指导下，编制程序在电子计算机上实际计算）1.利用现成的数学通用软件（如MATLAB、Mathematica、Maple等），对于N=32,64,128：(1）生成实数序列(x(k)-；(2）用FFT计算(x(k))-的离散Fourier变换序列(X(j))；(3)作出(x(k))和(IX()I)的图并进行分析（参见图16.5.4);(4）设定>0，将(IX(j))中满足/X()k%的数据全部置为零，再进行离散Fourier逆变换，将得到的数据与(x(k)比较；(5）改变s.的值，重复(4)，分析不同的s.对逆变换所得到的数据的影响。解源程序为function ex1601(N)t=0:N-1;x=randn (N, 1)*20;%randny=fft (x, N) ;z=abs(y);plot(t,x,"+",t,z,o')%delta=input（请输入误差）；for i=0:N-1if z(i+1)<deltay(i+1)=0;endendz=real(ifft(y));plot(t,x,'+",t,z,'o')运行结果分析：以N=128为例。本程序数据是随机产生的，“+”为原始数据，“o”为变换后的模的数据。4取=5，将X()中满足X()k的数据全部置为零，再进行离散6
计 算 实 习 题 （在教师的指导下，编制程序在电子计算机上实际计算） ⒈ 利用现成的数学通用软件（如 MATLAB、Mathematica、Maple 等）， 对于 N = 32, , 64 128 ： ⑴ 生成实数序列{ ( x k )}k N = − 0 1 ； ⑵ 用 FFT 计算{ ( x k )}k N = − 0 1 的离散 Fourier 变换序列{ ( X j)}N j= − 0 1； ⑶ 作出 { ( x k )}和{ | X j ( )| }的图并进行分析（参见图 16.5.4）； ⑷ 设定δ 0 > 0 ，将{ | X j ( )| }中满足 0 | X ( j)|< δ 的数据全部置为零， 再进行离散 Fourier 逆变换，将得到的数据与{ ( x k )}比较； ⑸ 改变δ 0 的值，重复⑷，分析不同的δ 0 对逆变换所得到的数据 的影响。 解 源程序为 function ex1601(N) t=0:N-1; x=randn(N,1)*20;%randn y=fft(x,N); z=abs(y); plot(t,x,'+',t,z,'o') % delta=input('请输入误差'); for i=0:N-1 if z(i+1)<delta y(i+1)=0; end end z=real(ifft(y)); plot(t,x,'+',t,z,'o') 运行结果分析：以 N=128 为例。本程序数据是随机产生的，“+”为 原始数据，“o”为变换后的模的数据。 取 0 δ = 5，将{ | X j ( )| }中满足 0 | X ( j)|< δ 的数据全部置为零，再进行离散 6
Fourier逆变换。“+”为原始数据，“o”为置零后变换得到的数据，与(x(k)比较几乎重合取。=50，同样处理后得到的数据，与(x(k))比较有些小误差。取8。=100，同样处理后得到的数据，与(x(k)比较误差清晰可见，但不很大。由于数据源不同，结果会有所差异。2.对于N=32,64,128，(1）产生两个实数序列(x(k))-和((k))-;(2)）用直接方法计算(x(k))和(y(k))的卷积(=(k)N=-;(3)改用离散Fourier变换的思想，用FFT计算(=(k))；(4)结合N比较两种算法所用的时间。解源程序为function t=ex1602(N)x=randn(N,1)*20;%randny=randn (N, 1)*20;%randntic%启动秒表%z=conv(x,y);
Fourier 逆变换。“+”为原始数据，“o”为置零后变换得到的数据， 与{ ( x k )}比较几乎重合 取 0 δ = 50，同样处理后得到的数据，与{ ( x k )}比较有些小误差。 取 0 δ =100，同样处理后得到的数据，与{ ( 比较误差清晰可见，但 不很大。 x k )} 由于数据源不同，结果会有所差异。 ⒉ 对于 N = 32, , 64 128 ， ⑴ 产生两个实数序列 { ( x k )}k N = − 0 1 和{ ( y k )}k N = − 0 1 ； ⑵ 用直接方法计算 { ( x k )}和{ ( y k )}的卷积{ (z k )}k N = − 0 1； ⑶ 改用离散 Fourier 变换的思想，用 FFT 计算{ (z k )} ； ⑷ 结合 N 比较两种算法所用的时间。 解 源程序为 function t=ex1602(N) x=randn(N,1)*20;%randn y=randn(N,1)*20;%randn tic %启动秒表 %z=conv(x,y); 7
for i=0:N-1z(i+1)=0;for j=o:iz(i+1)=z(i+1)+x(j+1)*y(i-j+1) ;endfor j=i+l:N-1z(i+1)=z(i+1)+x(j+1)*y (N+i-j+1);endendt1=toc;%计时tic;xl=fft(x,N) :yl=fft(y, N) ;zl=ifft(xl.*yl);t2=toc;t=[tl, t2];分析：计算所化时间与使用的计算机性能有关，由于计算机计时器的最小单位较大，对于较新的计算机，即使对于N=128，所化时间几乎为O。而且由于卷积采用代码解释执行速度较慢，Fourier变换采用内部函数速度很快，用FFT计算速度要快得多。武和的乘积，并与sin2x的3.用FFT计算多项式2n=o(2n+1)!=(2n)!Taylor级数的相应项比较。解源程序为function [z,maxerror]=ex1603(m);%z：乘积，maxerror：最大误差，m：阶数1en=4*m+2;a=zeros(len，1）；%被乘式系数a(2)=1;for i=4:2:2*m+2a(i)=-1*a(i-2) / (i-2) /(i-1) ;end8
for i=0:N-1 z(i+1)=0; for j=0:i z(i+1)=z(i+1)+x(j+1)*y(i-j+1); end for j=i+1:N-1 z(i+1)=z(i+1)+x(j+1)*y(N+i-j+1); end end t1=toc;%计时 tic; x1=fft(x,N); y1=fft(y,N); z1=ifft(x1.*y1); t2=toc; t=[t1,t2]; 分析： 计算所化时间与使用的计算机性能有关，由于计算机计时器的最 小单位较大，对于较新的计算机，即使对于 N=128，所化时间几乎为 0。而且由于卷积采用代码解释执行速度较慢，Fourier 变换采用内部 函数速度很快，用 FFT 计算速度要快得多。 ⒊ 用 FFT 计算多项式 ( ) ( )! − + + = ∑ 1 2 1 2 1 0 n n n m x n 和 ( ) ( )! − = ∑ 1 2 2 0 n n n m x n 的乘积，并与 sin 2 2 x 的 Taylor 级数的相应项比较。 解 源程序为 function [z,maxerror]=ex1603(m); % z:乘积，maxerror：最大误差，m：阶数 len=4*m+2; a=zeros(len,1);%被乘式系数 a(2)=1; for i=4:2:2*m+2 a(i)=-1*a(i-2)/(i-2)/(i-1); end 8
b=zeros（len，1）；%乘式系数b(1)=1;for i=3:2:2*m+1b(i)=-1*b(i-2)/(i-2)/(i-1) ;endc=zeros（len，1）%乘积系数c (2)=1;for i=4:2:lenc(i)=-4*c(i-2)/(i-1)/(i-2);endx=fft(a,len）；%Fourier变换y=fft（b，len）；%Fourier变换zl=x.*y;z=ifft（zl)；%Fourier逆变换maxerror=0;for i=l:lene=abs(z(i)-c(i));ife>maxerrormaxerror=e:endend9
b=zeros(len,1);%乘式系数 b(1)=1; for i=3:2:2*m+1 b(i)=-1*b(i-2)/(i-2)/(i-1); end c=zeros(len,1);%乘积系数 c(2)=1; for i=4:2:len c(i)=-4*c(i-2)/(i-1)/(i-2); end x=fft(a,len);%Fourier 变换 y=fft(b,len);%Fourier 变换 z1=x.*y; z=ifft(z1);%Fourier 逆变换 maxerror=0; for i=1:len e=abs(z(i)-c(i)); if e>maxerror maxerror=e; end end 9
计算结果误差分析：误差m10. 0520.00158732.756e-00543.006e-00752.248e-00961.224e-011随着m的增加，误差迅速减少。10
计算结果误差分析： m 误差 1 0.05 2 0.001587 3 2.756e-005 4 3.006e-007 5 2.248e-009 6 1.224e-011 随着 m 的增加，误差迅速减少。 10