习题164 Fourier变换和 Fourier积分 1.求下列定义在(-∞,+∞)的函数的 Fourier变换 A,00 (4)f(x)= x≥0, x=e a>0 0,x0)的正弦变换和余弦变换。 解正弦变换: f(o)= f(x)sin odx= e- sin oxdx= 余弦变换
习题 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
f(o)= f()cosoxdx= e cos oxdx +o 3.设f( rQJ(y=sinx,0≤rs3 xx≥0 求几*f2(x) 0,其它, 解记F()=f*(x)=*()1smO)0(x-0m,考虑∈D.2 当x≤0时,f(x-D=0,所以F(x)=0; 当x>z时,f(x-1)=,所以 F(x)=eesin(t)dt=e"(1+e2); 当0<x≤x时,f(x-0)= 所以 0, t F(x)= e- esin(Mhs、 sIn x-cosx+e 于是 ≤0 f *f2(x)=-(sinx-cosx +e),0<xs 2
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变换 1.说明离散 Fourier变换X(/)=∑x(n)e2可以看成 Fourier变换 f(o)」f(x)eatr 的离散近似形式的推广。 解假设 ∫f(x)eadx∑f(n)e 取Ax使 Ax 1 记W=,则k为整数时 于是 f()=∑W"∑f(N+n)△)x, 记x(n)=∑f(kN+n)Ax)Ax,所以 X()=f(0)=∑W"x(n)=∑x(n)e 2.证明正交关系式 解显然,j=k时,∑ 下面考虑j≠k,不妨设j<k。根据当≠1是方程xN=1的一个根时, 有∑n=0,令5 则ξN=e2 于是 k-y) 0。 3.设N=四(p,q∈N),构造只需O(p+q)N)次运算的 Fourier变换
习题 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
算法 解令W=eN,则k为整数时,W=W。假设 j=jq+j,方=0 l,j=0,1,…q-1, n=nP+n,n1=0,1,…,q-1,n=0,1,…,P-1。 X()=2xex以=∑x(mWm ∑xp+n)W W∑x(nP+n形 nJop 固定,计算xnp+mW需要q-1次乘法(n=0不需要做乘 法,对于相同的,∑x(mP+n1M是相同的,无需重复计算,所 有此类和式共需q(q-1)次乘法。对m求和需要(p-1)N次乘法,所以, 总共需要q(q-1)+(p-1)N=O(p+q)N)次乘法。 4.对N=23,具体写出以2为底的FFT的计算流程。 解记W=e8=e4,则W4=-1,W8=1。可得计算公式 X()=∑x(mWm,j=0,1…7 =[x(0)+(-1)x(4)+W[x(1)+(-1)x(5) +2/{[x(2)+(-1)yx(6)+W[x(3)+(-1)yx(7)]}。 计算流程 第一步: x1()=x(1)+x(+4) x1(+4)=W[x()-x(+4)]i=0,1,2,3
算法。 解 令 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
第二步: (1)=x()+x1(i+2) x2(i+2)=W2{x(i)-x(i+2)]i=0,1,45 第三步 X()=x2()+x2(i+1) X(+4)=x2()-x2(i+1,=0,2 X()=x2(+3)+x2(+4) X(+4)=x2(+3)-x2(+4),i=1,3
第二步： 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)}x; (2)用FFT计算{x(k)泌的离散 Fourier变换序列{X() (3)作出{x(k)}和{|X()}的图并进行分析(参见图16.5.4); (4)设定δ>0,将{XO)}中满足|X(kδ的数据全部置为零, 再进行离散 Fourier逆变换,将得到的数据与{x(k)比较 (5)改变δ的值,重复(4),分析不同的δ对逆变换所得到的数据 的影响。 解源程序为 function ex1601(N) t=0:N-1 x=randn(N, 1):*20; %rand y=fft(x, N) zp% lot(t,x,’+',t,z,’o') delta= input(请输入误差’) for i=0: N-1 f z(i+1)<delta y(i+1)=0 end z=real(ifft(y)) plot(t,x,’+’,t,z,’o’) 运行结果分析:以N=128为例。本程序数据是随机产生的,“+”为 原始数据,“o”为变换后的模的数据 °。889 取=5,将{X()}中满足|X(k的数据全部置为零,再进行离散
计 算 实 习 题 （在教师的指导下，编制程序在电子计算机上实际计算） ⒈ 利用现成的数学通用软件（如 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)}比较几乎重合 取a=50,同样处理后得到的数据,与{x(k)}比较有些小误差。 取δ=100,同样处理后得到的数据,与{x(k)}比较误差清晰可见,但 不很大 由于数据源不同,结果会有所差异 2.对于N=32,64,128 (1)产生两个实数序列{x(k)和{(k)} (2)用直接方法计算{x(k)}和{y(k)}的卷积{(k); (3)改用离散 Fourier变换的思想,用FFT计算{(k)} (4)结合N比较两种算法所用的时间。 解源程序为 function t=ex1602 ( N x=randn (N, 1)*20: %randn =randn ( N, 1)=20: %randn tic%启动秒表 %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 1=0: N-1 z(i+1)=0 for j=0: i z(i+1)=z(i+1)+x(j+1) end z(i+1)=z(i+1)+x(j+1)*y(N+i-j+1); end end t1=toc;%计时 xl=fft(x, N) yl=fft(y, N) zl=ifft(xl. *y1) t2=toc t=[t1,t2] 分析 计算所化时间与使用的计算机性能有关,由于计算机计时器的最 小单位较大,对于较新的计算机,即使对于N=128,所化时间几乎为 0。而且由于卷积采用代码解释执行速度较慢, Fourier变换采用内部 函数速度很快,用FFT计算速度要快得多 3.用FFT计算多项式∑x和∑x的乘积,并与2的 1)!m=0(2n Taylor级数的相应项比较 解源程序为 function [z, maxerror]=ex1603(m) %z:乘积, maxerror:最大误差,m:阶数 len=4*m+2 a= zeros(len,1);%被乘式系数 fori=4:2:2*m+2 a(i)=1*a(i-2)/(i-2)/(i-1) d
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; fori=3:2:2*m+1 b(i)=1*b(i-2)/(i2)/(i-1) 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变换 fft(b,len);% curler变换 z=ifft(z1);% Fourier逆变换 maxerror=O for i=l: len e=abs(z(i-c(i)) if e>maxerror maxerror-e d
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
计算结果误差分析: 误差 0.05 0.001587 3 2.756e-005 4 3.006e-007 2.248e-009 6 1.224e-01l 随着m的增加,误差迅速减少
计算结果误差分析： 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