数字信号处理(程佩青)课后习题解答(4)

第四章 快速傅立叶变换
运算需要多少时间。

计算需要多少时间，用，问直拉点的，用它来计算每次复加速度为平均每次复乘需如果一台通用计算机的FFT DFT[x (n)]512s 5 s 50.1μμ 解: 解: ⑴ 直接计算: 复乘所需时间:
复加所需时间: ⑵用FFT 计算:
复乘所需时间: 复加所需时间:
运算一次完成。

点试用一个为了提高运算效率值求今需要从值的点实序列是两个已知IFFT N n y n x k Y k X DFT n y n x N k Y k X ,,)(),()(),(,)(),()(),(.2s
N T N 01152.0 512log 105 log 105 2251262261=⨯⨯⨯=⨯⨯=--s T T T s
N N T 013824.0 002304.0 512log 512105.0 log 105.0 2126262=+=∴=⨯⨯⨯=⨯⨯⨯=--s
T T T s
N N T 441536.1 130816.0 )1512(512105.0 )1(105.0 21662=+=∴=-⨯⨯⨯=-⨯⨯⨯=--s N T 31072.1 512105 105 262
61=⨯⨯=⨯⨯=--
值的过程。

）（），（完成计算点）可用一次（）（）（综上所述，构造序列
）（）（）（）（可得：）（）（）（再根据都是实序列，
）（），（由原题可知：）
（）（）（）（（）（）（性质：
又根据可得序列点作对取序列依据题意解 ]Im[ ]Re[ ][][ ][ ).()( )()()( )()();()( :
:n y n x IFFT N k jY k X k Z n z n y n z n x n jy n x n z n y n x n jy n x k Y jIDFT k X IDFT k jY k X IDFT DFT n z IFFT N k Z k jY k X k Z k Y n y k X n x +===+=+=+=++=⇔⇔。

输出倒位序顺序频率抽取采用输入自然输出自然数顺序序时间抽取采用输入倒位流图抽取法的按时间抽取法及按频率画出基时), ,,( 2,16.3FFT N -=。

的运算量及乘不计乘相比较的基并就运算量与公式并画出流图基导出时)(2,4,16.4j j FFT FFT N ±±--=
)
,( )4()()()
,( )4()()(3
,2,1,03
,2,1,0,)( ,3
,2,1,03,2,1,0;, 440101011010102101011120102121k k X k k X k r k X k X n n x n n x n r n x n x k k k r k k N k k r r N n n n r n n N n r r N =+=+==+=+=∴⎩⎨
⎧==+=<=⎩⎨
⎧==+=<∴=⨯=有
对于频率变量同样令有对于解：依题意：
∑∑∑∑∑====++=+=
+==∴3
03
16
41641601
303
0)
4)(4(16
01
15
016
010
010******** )4( )4()()(n n k n k n k n n n k k n n n nk
W W W n n
x W n n x W n x k X
并画出流图。

采用基
的结果
算法求
为组合数时的
试用,)4
3
(
12
.5⨯
-
=
N
FFT
N
1 20 1
0 2
1
2
1
3,2,1,0
2,1,0 ,
,
,
4
3
r
r
N
n
n
n
r
n n
N
n
r r N
=
⎩
⎨
⎧
=
= +
=
<≤
∴
=
⨯
=
令
同样：有
对于
解：依题意：
∑∑∑
==
+
+
=
=
=
∴
=
+
=
+
=
=
+
=
+
=
∴
⎩
⎨
⎧
=
=
+
=
<
≤
3
2
)
3
)(
4(
12
1
11
12
1
1
1
1
1
1
2
1
1
1
1
01
1
1
)
,
(
)
(
)
(
)
,
(
)
3(
)
(
)
(
)
,
(
)
4(
)
(
)
(
2,1,0
3,2,1,0
,
)
0(
n n
k
k
n
n
n
nk
W
n
n
x
W
n
x
k
X
k
k
X
k
k
X
k
r
k
X
k
X
n
n
x
n
n
x
n
r
n
x
n
x
k
k
k
r
k
k
N
k
k有对于频率变量
并画出流图。

的结果同上题导出,52330.6⨯⨯==N
∑∑∑∑==+===∴=++==++=∴++===∴==
===
===
29
030
0120120120120123210301224
5010'
22103)3(30
0102010'
211
2
1
'1
010*********'102
3
1
2
0101)()()
,,()36()()
,,()510()()36( ),,( ),,()(4
,3,2,1,0,),,(),,( ),,(),,(1
,0,
),,(),,(),,(),,(2
,1,0,
),,(),,(02000111101202n nk
n k n n k k n k n k n n k n W n x k X k k k X k k k X k X n n n x n n n x n x k k k X k k k X k k k X k X k W n k k X
k k k X W n k k X n k k X k W
n n k X n k k X W n n k X n n k X k W
n n n x n n k X 则令
⎪⎩⎪
⎨⎧===++=++=<≤=⎪⎩⎪
⎨⎧===++=++=<∴=⨯⨯=2,1,01
,04
,3,2,1,0;36)300(,4,3,2,1,01
,02
,1,0 ;510 ,
5230
120120*********
120120313223
21k k k k k k k r k r r k k k k r r r N n n n n n n n r n r r n n N n r r r N 有对于频率变量令同样：有对于解：依题意：
流图如下图所示
:
M
N b M N a N z X DFT N N k e z N z
n x z X z n
n x n x n x M k N
j M n n
>≤-===⎩⎨
⎧≤≤=∙
∙∙∑-=- )(;
)(:
)(,
,1,,1,0,,)()( 0,1-M n 0 ),()(,
)(.721
之个抽样的方法，并证明的就能计算点找出用一个情况试对下列上的抽样即在上的抽样个等间隔点
在单位圆上变换我们希望计算求其他点的有限长序列研究一个长度为π
{}
2
00
01111010
2200
01020011
1210202101201205)3(3021063012204
030
330630530
15302
01
010300
1
2
4
2010)36)(510(30
1
2
4
0 ),,( ),,(),,(k n n k k k n n k n k n n n k n k n k n k n k n n n k n n n n k k k n n n n W W W W W n n n x W W W W W W
n n n x W
n n n x +========++++=⨯⎥⎥⎦
⎤⎢⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛=⨯=
=∑∑∑∑∑∑∑∑∑
∑∑
∑∑--=--=--=--=-+
+
+
=
<≤-=
≤∙
∙∙1
)1( 212 21
221
22)()()()()1()()( , )
(M N
l n k
n N
j N N n k
n N
j N n k
n N
j k N
j M n k
n N
j k
N
j e
n x e
n x e
n x e
X lN M N l e
n x e
X M N a ππππππ设依题意若解：
1
0,
)()(1
010)()(,)()()(, )( )]([)()1)1(0(]
)1([)()
10( ]
)2([)(
),()(),()(,
])1([)()(10
20
201
10
1
221210 2 )(21
)1(0
])1([21
)(21
2-≤≤=
⎩⎨
⎧-≤≤-≤≤=>≤=
∴---≤≤-+=-≤≤-+=+===-++
+
++
=
∑∑∑∑∑
∑∑-=--=-=-=----+----=-+--=+--=-∙
∙
∙
∙
∙
∙
∙
∙
∙
N k e n x
e
X N n M M n n x n x N n x M N b N z X DFT N n y
M N e
n y
e
X N l M n N l n x n y N n N l n x n y N n x n y n x n y e
e
e
N l n x e
N n x e
n x N n k
n N
j k N
j l m m
N n l m k
n N
j m
k N
j l l k
n N
j k
lN n N
j N l M n k
N l n j N n k
N n N
j N n k
n N
j πππππππππ则：即点到补零，可将若：点抽样
在单位圆上的即可计算点
，然后对它求一次，可先计算
由此可见，对于且令：
实现过程示意图。

的路径及画出平面路径为。

已知的复频谱点法求其前面
试用其他点序列已知一个CZT z A z z X CZT n x k k ; 20/2 ,2.1W ,3/ ,8.0)( 10 n
,07
n 0 ,1)( 8.80000πφπθ====⎩⎨⎧≤≤=
9
,,1,0,)]()([)(7
,,1,0 )(7,,1,0;)()(:)(}2.12.18.0{)2.1()(])([
90,)2.1()
8.0()()((1)
90,)
2.1(8.0 2.1;8.02
2 22
7
2
2
)(222
270
]
32)(20
220[2
)(2
1
22
220
22
2
222
1)310(70
7
)310(20
20300
=*======⨯⨯⨯=
∴--+=≤≤⨯=
=
∴≤≤⨯=====
=--=--
-=+-------+-=-=-+----∑∑
∑∑
k k h k g W
z X n
W
n h n W A n x n g W W
A n x W
e e z X n k k n nk k e z n x z X k e AW z e e
W W e e A A k k n n n n n k n n k n n n k n j n k n n k j k k n
k j n nk n
n n
k
k k j k
k
k j j j j 则：令则解：依题意：πππππ
ππππ
ϕπθ
由（1）式可得k z 的路径，如下表所示:
{}
时的抽样。

为实数
在变换不能计算即线性调频两者都不行两者都行
和为实数为实数使变换的的实轴上各点平面在点有限长序列计算一个可以用来变换线性调频的结论在下列说法中选择正确z H(z)z ,(b)(b)(a) (c)0
a ,a 1,-N ,1,
0,k ,z (b)1a ,a 1,-N ,1, 0,k ,)(:),( )()(..9k ≠==±≠==∙
∙∙∙
∙∙ak a z a z H z z z z n h M CZT z k k k
1,,1,0 , )()()()()()1( 0 ,0 , , 1 , )(,,,1,,1,0
)()()(102/2/)(102/2/10001000000)(001
0222200-=-===
±≠=====∴-====∙∙∙∙
∙∙∙∙∙
∑∑∑∑-=---=---=--+---=-N k n k p n g a a a n h a
a n h z H a a z a W A z z n h W A N k e W A AW z z
n h z H a M n k n k M n n
k M n nk k k
k k j k k k M n n k
k 此时即可
只需取
变换各点的上平面实轴
在若求有限长序列都是任意实数
其中：是正确的。

解：ϕθϕθϕθ
1
0,)()(,
)]([)( )
()()( ,
)]([)( ,
)]([)( 1
2 2/2-≤≤⋅==⋅===-+≥=-N k k r a z H L k R IFFT n r k p k G k R L n p FFT k p L n g FFT k G M N L FFT k k J 则点点点计算时可先求出
计算，式中取
为了用
10. 当实现按时间抽取快速傅立叶变换算法时，基本的蝶形计算
)
()()()()()(11q X W p X q X q X W p X p X m r N m m m r N m m -=+=++ 利用定点算术运算实现该蝶形计算时，通常假设所有数字都已按一定
比例因子化为小于1。

因此在蝶形计算的过程中还必须关心溢出问题。

(a) 证明如果我们要求2/1|)(|2/1|)(|<<q X p X m m 和
则在蝶形计算中不可能出现溢出，即
1
)](Im[, 1)](Re[1)](Im[, 1)](Re[1111<<<<++++q X q X p X p X m m m m
2/1|)](Im[| 2/1|)](Re[| 2
/1|)](Im[| 2/1|)](Re[| <<<<q X q X p X p X b m m m m ，，实际上要求）（
似乎更容易些，也更适合些。

问这些条件是否足以保证在蝶形 计算中不会出现溢出？请证明你的回答。

证明：(a)
2/1|)(| 2/1|)(|<<q X p X m m
1 |)(||||)(| |)()(||)(|1<⋅+≤+=+q X W p X q X W p X p X m r N m m r N m m 1
|)(||||)(| |
)()(||)(|1<⋅+≤-=+q X W p X q X W p X q X m r N m m r N m m
同理可证。

，证明故可用 )( )( 11q X p X m m ++
1|)(| 1|)(| 211<<++p X p X m m 故即
1)]([Im 1
)]([Re 1
)]([Im )]([Re 12121212<<<+++++p X p X p X p X m m m m 因此
1)](Im[1)](Re[11<<++p X p X m m 所以 1
)](Im[1
)](Re[11<<++q X q X m m 同理可证 )(b 证明：
2/1|)](Im[| 2
/1|)](Re[|<<p X p X m m 因为
2/1|)](Im[||)](Re[|4
/1|)](Im[|4
/1|)](Re[|2222<+<<p X p X p X p X m m m m 所以
2/2|)(|2/1|)(|2<<p X p X m m 或即。

不出现溢出证在蝶形运算中
所以上述条件不足以保。

不一定小于及
同理可得出。

不一定小于及的结果可以得出故利用，
不一定小于因此1|)](Im[||)](Re[:|1|)](Im[||
)](Re[|)(2/1|)(|1111q X q X p X p X a p X m m m m m ++++
()变换算法。

利用线性调频。

点快速傅里叶变换算法直接利用法的可能行：
的办法。

讨论下列各方然后再丢掉一些的抽样
用先算出比要求点数多个抽样。

计算时不能采的时在频率计算希望我们
，的傅里叶变换的有限长序列表示长度为z b a k k e X n x e X k j j )(10)( 10 )
9,,1,0(100/2)( )(10)(.112∙∙∙==πωωω
解：(a)
若直接利用10点快速傅立叶变换算法，则：
100290902)()()(n k j n n j n k j e
n x e n x e X k πωω-=-=∑∑==
将n 为偶数与n 为奇数的部分分开，可得：
9,,1,0, )()()12( )2()12(
)2()()()
( 11002050240100250
240100)12(2401002240
10021002222
22222∙∙∙=+=++=++=+
=
--=--=+-=-=--∑∑∑∑∑∑k k G e
k G e r x e
e r x e r x e
r x e n x e n x e X k j r k j r k j r k j r r k j r r
k j r n k j n n k j n k j ππππππππω为奇数为偶数 50250250240
222)2( )2( )2()(r k j r r k j r r k j r l e
l r x e
l r x e l r x k G πππ---=∑∑∑+++=+=
为奇数为偶数式中：
)
1,0( )24( )4( ))12(2(
))2(2(2521
050225220
)12(5021
0250220
22
222=++++=++++=-=--=+-=-=∑∑∑∑l e
l s x e e
l s x e l s x e
l s x s k j s k j s
k j s s k j s s k j s πππππ
(b) 如考虑利用线性调频z 变换算法，
则
∑∑=-=-==
9
0100/290))(( )()(n nk
k j n nk n k e n x W A n x z X π 变换算法。

所以不能利用线性调频的函数是的函数。

因为这里必须不是在应用这种算法时，z k W k
W ,
12. 我们希望利用一个单位抽样响应为N=50个抽样的有限冲激响应滤波 器来过滤一串很长的数据。

要求利用重叠保留法通过快速傅立叶变换 来实现这种滤波器，为了做到这一点 ,则:
(1) 输入各段必须重叠P 个抽样点 ；
(2) 我们必须从每一段产生的输出中取出Q 个抽样点，使这些从每一段得
到的抽样连接在一起时，得到的序列就是所要求的滤波输出。

假设输 入的各段长度为100个抽样点，而离散傅立叶变换的长度为128点。

进一步假设，圆周卷积的输出序列标号是从n=0到n=127。

则：(a)求P ； (b)求Q ； (c)求取出来的Q 个点之起点和终点的标 号，即确定从圆周卷积的128点中要取出哪些点，去和前一段的 点衔接起来。

解：
(a) 由于用重叠保留法，如果冲激响应
h(n) 的点数为N 点，则圆周卷积
结果的前面的(N-1)个点不代表线
性卷积结果。

故每段重叠点数P 为
P=N – 1 =50 – 1=49
(b) 每段点数为 27 =128，但其中只
有100个是有效输入数据，其余
28个点为补充的零值点。

因而
各段的重叠而又有效的点数Q为
Q=100 – P=100 – 49 =51
(c) 每段128 个数据点中，取出来的
Q个点的序号为n=49 到n=99。

用这些点和前后段取出的相应点
连接起来，即可得到原来的长输
入序列。

另外，对于第一段数
据不存在前一段问题，故在数据
之前必须加上P=N – 1 =49个
零值点，以免丢失数据。

13. 请用C语言编写程序:
(1) 按频率抽取的FFT算法(2) 分裂基FFT算法解：（1）
/*Free_Copy*/
/* C语言编写的频率抽取FFT算法(最大计算64点）*/
/* 输入: 序列点数、序列值* /
/* 输出: 序列FFT变换后的数值及反变换(应与原序列相同) */
#include "conio.h"
#include "math.h"
#include "stdio.h"
#define N 64
#define PI 3.1415926
#define w0 (0.125*PI)
#define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x;
#define Cequal(a,b) a.x=b.x;a.y=b.y;
#define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y;
#define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y;
#define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n);
struct comp
{
float x;
float y;
};
void main()
{
int i,j,nu2,nm1,n,m,le,le1,k,ip,z;
int flag,f,n1;
struct comp a[N],t,t1,w,d;
float a_ipx,m1;
printf("\nThis program is about FFT by DIF way. ");
printf("\nplease enter N : ");
scanf("%d",&n1);
n=n1;
m1=log(n1)/log(2);
m=log(n1)/log(2);
if (m!=m1) n=pow(2,m+1);
for(i=0;i<n;i++) {a[i].x=a[i].y=0.0;}
printf("\n");
for(i=0;i<n1;i++)
{
printf("\nplease enter data(%d)_[Re]: ",i); scanf("%f",&a[i].x);
printf("\nplease enter data(%d)_[Im]: ",i); scanf("%f",&a[i].y);
}
for(z=0;z<=1;z++)
{
flag=-1;
for (m=(log(n)/log(2));m>=1;m--)
{
le=pow(2,m);
flag++;
le1=le/2;
for( j=0;j<le1;j++)
{
for (i=j;i<=(n-1);i+=le)
{
ip=i+le1;
Cequal(t,a[i]);
Cequal(t1,a[ip]);
f=(int) (i*pow(2,flag))%n;
Wn(w,f);
Cadd(a[i],t,t1);
Csub(a[ip],t,t1);
a_ipx=a[ip].x;
if (z==1)
{
w.y*=-1;
}
a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
a[ip].y=a_ipx*w.y+a[ip].y*w.x;
}
}
}
nu2=n/2;
nm1=n-2;
j=0;i=0;
while(i<=nm1)
{
if (i<j)
{
Cequal(d,a[j]);
Cequal(a[j],a[i]);
Cequal(a[i],d);
}
k=nu2;
while(k<=j)
{
j=j-k;k=k/2;
}
j=j+k;
i=i+1;
}
if(z==0)
{
printf("\n序列的fft是:\n\n");
}
else
printf("\n用ifft计算出的原序列是:\n\n" ) ; for(i=0;i<n;i++)
if(z==0)
{
printf(" %7.3f",a[i].x);
if (a[i].y>=0)
printf(" + %7.3f j \n",a[i].y);
else
printf(" - %7.3f j \n",fabs(a[i].y));
a[i].y= -a[i].y;
}
else
{
printf(" %7.3f",a[i].x/n);
a[i].y=-a[i].y/n;
if (a[i].y>=0)
printf(" + %7.3f j \n",a[i].y);
else
printf(" - %7.3f j \n",fabs(a[i].y));
}
}
printf("\n");
}
（2）;分裂基FFT 算法程序
/*Free_Copy*/
/*主程序:64点分裂基FFT算法*/
/*输入: 64点任意序列*/
/*输出: 序列的FFT变换*/
#include "conio.h";
#include"math.h"
#include"stdio.h"
#define PI 3.1415926
#define N 128
void main()
{
float x[N],y[N],xt;
float cc1,cc3,ss1,ss3;
float r1,r2,r3,s1,s2,a,a3,e,m1;
int n,n1,m,j,k,i;
int is,id,i0,i1,i2,i3,n2,n4;
printf("\nThis program is about FFT by SPEFT way. ");
printf("\nplease enter n : ");
scanf("%d",&n1);
n=n1;
m1=log(n1)/log(2);
m=log(n1)/log(2);
if (m!=m1) n=pow(2,m+1);
for(i=0;i<=N;i++)
{
x[i]=y[i]=0.0;
}
printf("\n");
for(i=1;i<=n1;i++)
{
printf("\nplease enter data(%d)_[Re]: ",i);
scanf("%f",&x[i]);
printf("\nplease enter data(%d)_[Im]: ",i);
scanf("%f",&y[i]);
}
j=1;
for (i=1;i<=n-1;i++)
{
if (i<j)
{
xt=x[j];
x[j]=x[i];
x[i]=xt;
xt=y[j];
y[j]=y[i];
y[i]=xt;
}
k=n/2;
while (k<j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
is=1;
id=4;
while (is<n)
{
for (i0=is;i0<=n;i0+=id)
{
i1=i0+1;
r1=x[i0];
x[i0]=r1+x[i1];
x[i1]=r1-x[i1];
r1=y[i0];
y[i0]=r1+y[i1];
y[i1]=r1-y[i1];
}
is=2*id-1;
id=4*id;
}
n2=2;
for (k=2;k<=m;k++)
{
n2=n2*2;
n4=n2/4;
e=2.0*PI/n2;
a=0.0;
for (j=1;j<=n4;j++)
{
a3=3.0*a;
cc1=cos(a);
ss1=sin(a);
cc3=cos(a3);
ss3=sin(a3);
a=j*e;
is=j;
id=2*n2;
while (is<n)
{
for (i0=is;i0<=n-1;i0+=id) {
i1=i0+n4;
i2=i1+n4;
i3=i2+n4;
r1=x[i2]*cc1+y[i2]*ss1;
s1=y[i2]*cc1-x[i2]*ss1;
r2=x[i3]*cc3+y[i3]*ss3;
s2=y [i3]*cc3-x[i3]*ss3;
r3=r1+r2;
r2=r1-r2;
r1=s1+s2;
s2=s1-s2;
x[i2]=x[i0]-r3;
x[i0]=x[i0]+r3;
x[i3]=x[i1]-s2;
x[i1]=x[i1]+s2;
y[i2]=y[i0]-r1;
y[i0]=y[i0]+r1;
y[i3]=y[i1]+r2;
y[i1]=y[i1]-r2;
}
is=2*id-n2+j;
id=4*id;
}
}
}
printf("\n分裂基fft结果是：\n ");
for (i=1;i<=n;i++)
{
printf("\n %7.3f, %7.3fj",x[i],y[i]);
y[i]=-y[i];
}
getch(); printf("\n\n"); }。