数值分析 - 第5章 快速傅里叶变换

F (ω ) =
1 α 2 +ω2 ω α
ϕ (ω ) = −arctg (
ω
4．解： X (ω ) = 1 2π
∫
∞
−∞
exp( −iω t ) x (t ) dt =
1 2π
∫
T1
−T1
exp( −iω t ) dt =
sin ωT1 πω
为使谱分辨率 提高一倍，F=5Hz，要求
5.4 习题解答或提示 1 略。 2 求对称指数函数 f(t)的傅立叶变换 f (t ) = e
=
1 2π
∫
T
−T
exp(−iω t )hdt =
代入傅立叶积分公式，得
F[ F (t )] = F [
f (t ) = ∫
∞
−∞
h sin ω T exp(iω t ) dω πω
2π ω < 1 2 sin t ] = 2πf (ω ) = t 0 ω >1
例 2 已知 x(4)=R4(n), 分别求 N=8 和 N=16 时的 X(k) 解： N＝ 8 时
图 5－3 例 7 令 x(n)=(-0.9) , n=[-5:5],讨论其离散时间傅立叶变换的 共轭对称性。 解：n=[-5:5]; x=(-0.9).^n; k=-200:200; w=(pi/100)*k; X=x*(exp(-j*pi/100)).^(n'*k); magx=abs(X); angx=angle(X); subplot(2,1,1); plot(w/pi,magx); grid; title('幅度 部分') subplot(2,1,2); plot(w/pi,angx/pi); grid; title('相角部分') 如图 5－4 所示
解:
如图 5－7 所示
Fra Baidu bibliotek亦即
|X(Ω)|= |X(－Ω)|，φ (Ω)= －φ( －Ω)
所以描述实信号频谱分布时 只要考虑
半 个周期数字频率范围
例 14 利用 DFT 的矩阵表达式求 4 点序列 x (n ) = c o s ( 2π n/ N )的离散傅立叶变换
解:由 N=4,得
于是
3．解：
F (ω ) = ∫
∞
−∞
f (t )e − jω t dt =
1 α + jω
(α > 0)
幅频 相频 例 15 对 实信号进行谱 分析 ，要 求 谱分 辨率 F≤ 10Hz， 信号 最高 频率 fc=2.5KHz,试确定 最小纪录时间 Tpmin，最大的采样间隔 Tmax，最少 的采样点数 Nmin。如 果 fc 不变，要 求谱分辨率增加一 倍，最少的 采样点数和最小的纪录时间是多少？ 解: 因此， Tpmin=0.1s 。因为要求 fs ≥2fc， 所以
n
图 5 －4 例 8 比较序列 x(t)=cos(80 π t)分别经过 fft 和 ifft 变换后的序列 解： MATLAB 命令窗口输 入如下 n=0:127;N=128 t=0:0.01:pi; %时域数据点 q=n*2*pi/N; %频域 数据点数 X=cos(80*pi*t); Y=real(ifft(fft(X))); subplot(2,1,1); plot(t,X) title(' 原信号') subplot(2,1,2); plot(t,Y) title('两次 变换后恢复的 信号') 如图 5－5 所示
π 4k − j2 16 π − j2 k 16
1− e
=e
πk − j3 16
sin(πk / 4) sin(πk / 16)
（k=0,1,2 ….15）
由该例可知：频率采样点数不同， DFT 的长度不同，DFT 的结果也不同。 例 3 试求函数 sin( t ) 的傅立叶变换 t
解：若直接利用公式
来求傅立叶变换将不是一件容易的事情。这里利用对称法求解。
f (t ) ↔ Eτ
sin
1 f (t ) = 0
t ≤1 t 1
图 5-1 例 5 给定数 学函数 x(t)=12sin(2π×10t+π/4)+5cos(2π×40t) 取 N=128，试对 t 从 0~1 秒 采样，用 fft 作快速傅立叶变换，绘制相应的振幅 -频率图。 在 0~1 秒时间范围内采样 128 点，从而可以确定采样周期和采样频率。由于离散傅立叶变换时的下标应 是从 0 到 N-1 ，故在实际应 用时下标应该 前移 1。又考虑到对离散傅立叶变换来说，其振幅 | F(k)|是 关 于 N/2 对称的，故只须使 k 从 0 到 N/2 即可 程序如下： N=128; % 采样点数 T=1; % 采样时间终点 t=linspace(0,T,N); % 给出 N 个采样时间 ti(I=1:N) x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t); % 求各采样点样本值 x dt=t(2)-t(1); % 采样周期 f=1/dt; % 采样频率 (Hz) X=fft(x); % 计算 x 的 快速傅立叶变换 X F=X(1:N/2+1); % F(k)=X(k)(k=1:N/2+1) f=f*(0:N/2)/N; % 使频率轴 f 从零 开始 plot(f,abs(F),'-*') % 绘制振幅- 频率图 xlabel('Frequency'); ylabel('|F(k)|') 如图 5－2 所示
图 5－2 例 6 对一给定 的连续信号 2sin(4 πt )+5cos(8 πt )画出其对应 的频谱图形 (N=64) 解：N=64; n=0:N-1; t=0.01*n; q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=64') 如图 5-3 所示
=
N=16 时
1− e
π − j 28 4k πk − j 28
1− e
=e
πk − j 38
sin(πk / 2) sin(πk / 8)
2πkn −j 16
X (k ) = ∑ x (n)WNnk = ∑ R4 (n)e
n =0 n =0
N −1
15
= ∑e
n=0
3
π kn − j2 16
=
1− e
−α t
F (ω ) =
2α α 2 +ω2
5．编写一个 M 文件，代码 如下： function y=dft(x,N) %x 为 输入信号 %N 为所取的点数 y=zeros(1,N); for k=1:N for m=1:N y(k)=y(k)+x(m)*exp(-j*2*pi*(k-1)*(m-1)/N); end; end;
（a） ；在命令窗口 中运行 x=[1 1 -1 -1],dft(x,4); 结果如下： 0 2.0000 - 2.0000i 0 2.0000 + 2.0000i （b） ：在命令窗口中运 行 y=[1 2 3 4 ]; dft(y,4); 结果如下： 10.0000 -2.0000 + 2.0000i -2.0000 - 0.0000i -2.0000 - 2.0000i 6．提示：参 考例 6 的解答 7．提示：参 考例 10 的解答 8．提示：参考例 10 的解答 9．解 在命令窗口 中输入 x=[1 2 3 4 0 1 2 3] 结果如下： Columns 1 through 4 16.0000 1.0000 - 2.4142i -4.0000 + 4.0000i 1.0000 - 0.4142i
每一周期取样点数 N=T0 / T
s1
=16，则
(b) 同理可求得当 fs=8 样点/周期时
小结 : （ 1）离散时间周期信号的频谱 X ( kΩ0 ) 是具有谐波性的 周期序列，而连续时间周期信号的频谱 X(k ω0)是 具有谐波性的 非周期序列。X( kΩ0 ) 可以 看作 X(kω0)的 近似式，近似 程度与取样间隔 T 值 的选取 有关 。 （ 2）在满足 取样定理的条 件下，从一个 连续时间频带有限的周期信号得到的周期序列，其频谱在 | Ω|<π 或 | f | < (f s / 2 ) 范围内等于原始信号的离散频谱。因此可以利用数值计算的 方 法，通 过计算机，方便的截取一个 周期的样点 x(n) ，准确的求出连续信号 的各谐波分量 X（k ω0）。 （3）在不满足 取样定理条件 下，由于 X ( kΩ0 ) 出现频谱 混叠，存 在混叠误差。当误差比较小时 ，X ( k Ω0 ) 在 0< f < ( f s / 2 )范围内可作为 X（kω 0） 的近似；当误差较大， 则无法从 X ( kΩ0 )求 得 X（kω 0） 。为此，必须减小取样间隔， 增加取样率， 尽量使 X ( kΩ 0 )逼近 X（kω0 ） 。 例 11 求 非周期序列 x (n )= a ε(n )，a=0.8 的傅立叶变换。 解： 解: 幅度频谱与相 位频谱分别为
5．3 习题课 例题 1 将矩形脉冲 f (t) = h rect(t/2T)展开为傅立叶积分。 解：先求出 f (t) 的傅立叶变换
F (ω ) =
↔ F (ω ) = F[ f (t )] =
h sin ωT πω
根据偶函数对称原理得
2 sin ω ω
∫ 2π
1
∞
−∞
exp( − iω
t ) f ( t ) dt
n
图 5－7 例 13 已知一周期连续频谱如图 5－8 所示，求其相应的序列 。
结论
从以上三个例子 可以总结出，离散时间非周期实信号 的频谱是连续 的周期性频谱，具有埃尔米特 例 12 求 有限长序列 x (n ) 的频谱并作图 ，已知 x(n)=1( —M≤n≤M ),M=2,其余为 0 (Hermitian symmetry) 对称性质， 即幅度对频率 具有偶对称，相位具有奇对称，表示为
F (ω ) = ∫ f (t )e − jωt dt
−∞
∞
=∫
ωτ 2 ωτ 2
sin t − jωt e dt −∞ t
∞
w1=[0:1:500]; w=(pi/500)*w1; x=exp(i*w)./(exp(i* w)-0.5*ones(1,501)) magx=abs(x); angx=angle(x); realx=real(x); imagx=imag(x); subplot(2,2,1);plot(w/pi,magx); xlabel('以 pi 为单位的频率 ');title('幅度部分 ');ylabel('幅度 ') subplot(2,2,3);plot(w/pi,angx); xlabel('以 pi 为单位的频率 ');title('相角部分 ');ylabel('弧度 ') subplot(2,2,2);plot(w/pi,realx); xlabel('以 pi 为单位的频率 ');title('实部 ');ylabel('实部') subplot(2,2,4);plot(w/pi,imagx); xlabel('以 pi 为单位的频率 ');title('虚 ');ylabel('虚部') 如图 5-1 所示
= 2 sin( 6πt ) + 7 cos( 8πt ) 求 N=128 点 DFT 的幅度谱和相位谱。
解：MATLAB 命令窗口输入如下： N=128; n=0:127; t=0:0.01:127; %时 域数据点数 q=n*2*pi/N; %频 域数据点数 X=2*sin(6*pi*t)+7*cos(8*pi*t); Y=fft(X,N); subplot(2,1,1); plot(q,abs(Y)) subplot(2,1,2); plot(q,angle(Y)) 如图 5－6 所示
图 5－5 例 9 已知模拟信号 x (t )
图 5－6 例 10 已知一连续时间 周期信号 x ( t )=2 cos 6πt + 4 sin 10πt ,现以不同的取样率 (a) fs1 =16 样点 /周期 (b ) f s2=8 样点/周期， 对它进行取样。 试分别求出取样后周期序列的频 谱并于原始信号的频谱作一 比较。 解:(a) 按 fs1=16,Ts1=1/16 所以
F[
3 − j 28π kn
sin t π ω < 1 ]= t 0 ω > 1
X (k ) = ∑ x(n)WNnk = ∑ R4 (n)e
n =0 n =0
N −1
7
−j
2 πkn 8
= ∑e
n =0
= π [u (ω + 1) − u (ω − 1)
例 4 将 [0,π ] 分为 501 个等间隔的点，计算 x(n)=(0.5)nu(n) 的离散傅立叶变换 X ( e jϖ ) ，并画出其模、相角 、实部、虚部的 曲线。 解：MATLAB 程序如下：