Chơng
6
Biến đổi Fourier rời rạc
6.1 Chỉ dẫn
Trong chơng 2,chúng ta đã chứng minh rằng đáp ứng tần số của hệ thống của hệ thống
tuyến tính bất biến (LSI ) 2-D đợc cho bởi:


=

=
+
=
1 2
2211
)(
2121
),(),(
k k
kkj
ekkhH


(6.1)
Nếu h(k
1
,k
2
) chỉ có chỉ tồn tại với k
1




(6.2)
Công thức này chứng tỏ rằng
H( , )

1 2
là tuần hoàn, chu kỳ tuần hoàn là 2. Nếu
chúng ta lấy mẫu dới dạng
1
,
2
, và miền xác định là (0
1
2) và (0
2
2), N
ì
N mẫu, chúng ta có thể viết:



1 1
2
=
N
n
và



N
j
N
k
ekkhnnH

(6.4)
Biểu thức (6.4) đợc gọi là biến đổi Fourier rời rạc 2-D hay còn gọi là DFT. Công thức
này đợc áp dụng vào nhiều ứng dụng nh lọc, nén ảnh, phóng đại ảnh. Trong chơng này
chúng ta sẽ nghiên cứu 2-D DFT và các kỹ thuật tính toán. Đầu tiên, chúng ta sẽ xem xét
1-D DFT, sau đó mở rộng ra cho 2-D.
6.2 Biến đổi Fourier 1-D
Biến đổi Fourier 1-D cho tín hiệu thời gian rời rạc f(kT) tính theo công thức :



=

=
1
0
2
)()(
N
k
nk
N
j
ekTfnF


)()(
nj
enAnF

=
(6.7)
Ký hiệu A(n),

(n) gọi là phổ khuyếch đại và phổ pha của F(n).
6.2.1 Biến đổi ngợc DFT
Hàm f(k) là biến đổi ngợc DFT của F(n) cho bởi theo biểu thức



=
=
1
0
2
)(
1
)(
N
n
nk
N
j
enF
N
kf

0
)(
1
0
1
0
1
0
)(
1
)(
1
)(
1
N
m
N
n
mkn
N
kn
N
N
n
N
m
nm
N
N
n

(k -m )
+ W
N

2(k -m )
+ + W
N

(N-1)(k -m )

hoặc

)(
2
))(2(
m)-(k
N
m)-N(k
N
1
1
W-1
W-1
mk
N
j
mkj
e
e
S

N
==


=
Kết quả này giống nh biểu thức (6.8).
Khi f(k) có thể rút ra từ F(n) và ngợc lại, chúng gọi là cặp biến đổi. Cặp biến đổi này
có dạng

)()( nFkf
76
Chú ý từ biểu thức (6.8) ta có thể dễ dàng chứng minh:

)(
)(
1
)(
1
1
0
.
2
1
0
)(
2
kf
enF
N
enF

1
(n) và f
2
(n), và cả hai dãy này tuần
hoàn với chu kỳ N, đợc dùng để tính
f
3
(k) = af
1
(k) + bf
2
(k) (6.11)
là kết quả của biến đổi DFT f
3
(n) cho bởi
F3(n) = aF1(n) + bF2(n) (6.12)
ở đây a, b là hằng số và
F
1
(n) = DFT của f
1
(k)
F2(n) = DFT của f2(k)
Tính đối xứng. Tính đối xứng của DFT rất hay đợc dùng.
nk
N
j
N
k
nk



2
1
0
2
1
0
2
1
0
)(
)(
)(
)()(
Nếu f(k) là thực thì




=

=








2
)
77
)( Nkf +
khi
)()(
1
0
1111
1
11


=

=
N
k
kn
N
WkfnF

)()(
1
0
2222
2
22




=
=
1
0
213
)().(
1
)(
vì vậy









=








=


)(
2211
21
21
1 2
21
1
)()(
)()(
1
N
n
kkkn
N
N
k
N
k
nk
N
N
n
N
k
N
k
kkn
N
3
W

01
113
lNkkfkfkf
N
k
=


=
(6.14)
ở đây k = 0 đến 2N - 1.
Biểu thức trên biểu diễn tích chập của hai tín hiệu tuần hoàn. Chú ý rằng biêủ thức này
chỉ áp dụng cho hai dãy có chung một chu kỳ, và chiều dài của dãy tính theo biểu thức
trên là 2N - 1. Kết quả này chứng minh rằng trong DFT, tín hiệu có số mẫu lớn hơn N sẽ
đợc biến đổi thành dãy tuần hoàn có chu kỳ N. Khi dùng DFT cho một tín hiệu không có
chu kỳ, mà kết quả thu đợc từ tích hai dãy, ta sẽ phạm một sai lầm gọi là lỗi wraparound.
Đó là lý do ta phải làm cho cả hai dãy có chu kỳ bằng nhau. Để sửa lỗi này, một số số 0
cần phải thêm vào cả hai dãy để chiều dài hai dãy bằng nhau. Ví dụ, nếu một dãy có
chiều dài A, một dãy có chiều dài B, kết quả ta phải thêm các số 0 cho cả hai dãy có
chiều dài ít nhất là A + B - 1.
78
cho k = k
1
+ k
2
+ lN
các trờng hợp còn lại.
W =
.=
2

N
k
kn
N
N
k
kn
N
Wkfkf
WkfWkfnFnF



=


=

=


=

Bài tập 6.1
Cho hai dãy sau




=

2
(k) để chu kỳ của chúng = 5 + 6 - 1. Làm
lại phần 3 và so sánh kết quả.
6.3 Thuật toán biến đổi nhanh Fourier
Tính trực tiếp giá trị của DFT bao gồm N phép nhân phức và N - 1 phép cộng phức
cho mỗi giá trị của F(n). Khi N giá trị đợc tính toán thì N
2
phép nhân và N(N - 1) phép
cộng đợc tính toán. Cũng nh vậy, cho N có giá trị rất lớn, tính trực tiếp giá trị của DFT sẽ
đòi hỏi một số phép tính lớn đến mức không thể chấp nhận đợc. Để ví dụ, cho N = 1024
= 2
10
ta sẽ phải tính 2
20
= 1,048,576 phép nhân số phức và một số gần bằng nh vậy các
phép cộng.
Hoàn thiện có nghĩa là phải giảm số phép tính trong biến đổi Fourier xuống. Dới đây
chúng ta sẽ giới thiệu hai thuật toán hay dùng là thuật toán phân chia thời gian và thuật
toán phân chia tần số. DFT dùng các thuật toán trên gọi là Fast Fourier transform (FFT).
6.3.1 Thuật toán phân chia thời gian
Xem xét tính toán của DFT cho bởi (5.6) với N= 2
r
(r là một số nguyên bất kỳ). Cơ sở
của thuật toán phân chia thời gian thì rõ ràng. Tuy nhiên, việc thiết kế phần mềm cũng
đòi hỏi một số phân tích chi tiết. Để làm rõ các bớc của thuật toán này chúng ta sẽ bắt
đầu phân tích với N = 16 và sau đó mở rộng ra áp dụng cho N bất kỳ.
Cơ sở của thuật toán phân chia thời gian dựa trên cơ sở chiến lợc chia và chiếm. Các
bớc sau sẽ làm sáng tỏ thuật toán. Vì trong trờng hợp này N =16; nên,



nk
WkfWkfnF

k chẵn k lẻ
Chúng có thể viết thành
79
0 k
1
4
các trờng hợp còn lại.
0 k
1
5
các trờng hợp còn lại.

=
+
=

++=
7
0
)12(
16
7
0
)2(
16
)12()2()(
k


=

++=
7
0
816
7
0
8
)12()2()(
k
nkn
k
nk
WkfWWkfnF
đặt

f(2k)(k)f
10
=

1)f(2k(k)f
11
+=
Ta đợc

=

=


=
7
0
81111
)()(
k
nk
WkfnF
(6.18)
Viết lại biểu thức (6.16) chúng ta đợc

(n)FW (n)F F(n)
11
-n
1610
+=
(6.19)
Cũng nh vậy, phát triển cho một biểu thức
(n)FW -(n)F 8)F(n
11
-n
1610
=+
(6.20)
Biểu thức (6.19) và (6.20) định dạng những đơn vị tính toán gọi là bớm. Hình 6.1 là
biểu đồ của phần tử bớm. Ký hiệu W
16
-n
thờng gọi là trọng lợng hay hệ số xoay. Hai biểu

k
nkn
k
nk
WkfWWkfnF
DÔ thÊy

-2n
16
-n
8
WW =
®Æt

1)(2kf(k)f
(2k)f(k)f
1021
1020
+=
=

H×nh 6.1 (a) Bím; (b) BiÓu diÔn rót gän.
81
F
10
(n) F(n)

F
11
(n) F(n+8)

nk
k
nnk
WkfWWkfnF

(n)FW (n)F (n)F
21
-2n
162010
+=
(6.21)

(n)FW - (n)F 4)(nF
21
-2n
162010
=+
(6.22)
ở đây

=

=
3
0
42020
)()(
k
nk
WkfnF


=
3
0
42222
)()(
k
nk
WkfnF
(6.27)
82
012!4567
012!4567
012!4567
891011121!
1415
F
10
(n)
F(n)
F
11
(n)
X(k)
X(k)

=

=
3

(6.30)

(n)FW (n)F (n)F
33
-4n
163221
+=
(6.31)
(n)FW - (n)F 2)(nF
33
-4n
163221
=+
(6.32)
(n)FW (n)F (n)F
35
-4n
163422
+=
(6.33)
(n)FW - (n)F 2)(nF
35
-4n
163422
=+
(6.34)

(n)FW (n)F (n)F
37
-4n


(6.38)
, vv.
Các biểu thức từ (6.29) đến (6.36) cho kết quả trong bớc thứ ba của thuật toán và biểu
diễn trong lu đồ hình 6.4.Mỗi phần tử từ F
30
(n) đến F
37
(n) có thể chia tiếp thành hai phần
tử nữa và bớc này tạo thành sơ đồ cuối cùng (bớc đầu tiên) trong lu đồ.
8!
012!012!
012!012!
012!4567
012!4567
0
F
20
(n)
2
F
21
(n)
4
6
F
22
(n)
0
F

phép nhân phức là 8/2 . 4. Tổng quát cho N = 2
r
số phép nhân phức là (N/2) . r = (N/2 )
log
2
N và số phép cộng là Nlog
2
N. Chú ý, thực tế số phép nhân sẽ giảm xuống một ít, vì
trong bớc đầu tiên hệ số xoay W
0
= 1 và trong các bớc còn lại chúng ta cũng có các bớm
với hệ số xoay = 1.
Xem xét trờng hợp N = 1024 = 2
10
. Số phép nhân cần dùng cho FFT là (N/2).10 =
1024 ì 5 = 5120 so với 1 triệu phép nhân cho tính trực tiếp biến đổi DFT, đây là phơng
pháp tiết kiệm thực sự cho tính toán.
Bây giờ, chúng ta sẽ vạch ra thuật toán FFT. Đó không đơn thuần chỉ là sự phát triển
một chơng trình từ lu đồ. Tuy nhiên, chúng ta có thể nghiên cứu lu đồ và vạch ra các bớc
có thể dùng để phát triển một chơng trình. Từ lu đồ của hình 6.5 chúng ta có thể viết:
Bớc thứ nhất. Trong bớc này ta có tám bớm với trọng lợng (hệ số xoay) W
0
= 1.
Chúng ta có thể viết (xem hình 6.6)

for (j=0 đến 15 với bớc tăng 2)
{
T=X(j+1);
X(j+1)=X(j) - T;
X(j)=X(j) + T;

F
!6
(n)
F
!7
(n)
X(k)
X(k)
Dãy
đầu
vào
đã
đ ợc
sắp
xếp
lại
Bớc thứ hai. Chúng ta có:
1.Bốn bớm với trọng lợng bằng 1.
for (j=0 đến 15 với bớc tăng 4)
{
T=X(j);
X(j+2)=X(j) - T;
X(j)=X(j) + T;
}
2. Bốn bớm với trọng lợng W
4
= W(3). Chú ý rằng chúng ta coi rằng các hệ số xoay
W, W
2
, , W

3. Hai bớm với trọng lợng bằng W(3) = W
4
.
for (j=2 đến 15 với bớc tăng 8)
{
87
0
1
2
3
4
5
6
7
0
2
4
6
0
2
4
6
0
4
0
4
0
4
0
4

0110
1110
0001
1001
0101
1101
0011
1011
0111
1111
0
8
4
12
2
10
6
14
1
9
5
13
3
11
7
15
0
4
8
12

2
3
4
5
6
7
8
9
10
11
12
13
14
15
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

.
for (j=3 đến 15 với bớc tăng 8)
{
T=X(j)W(5);
X(j+4)=X(j) - T;
X(j)=X(j) + T;
}
Bớc thứ t và là bớc cuối cùng.
1.Một bớm với trọng lợng bằng 1.
T = X(0)
X(8)= X(0) - T
X(0) = X(8) +T
2. Một bớm với trọng lợng bằng W

(0) = W.T = X(1)W(0)
X(1+8)= X(1) - T
X(1) = X(1) +T
3. Một bớm với trọng lợng bằng W

(1) = W
2
.T = X(1)W(1)
X(2+8)= X(0) - T
X(2) = X(2) +T

Các bớc dẫn chúng ta đến thuật toán với N = 16.
Thuật toán
ip=1
kk=8
incr=2
cho iter=0 đến 3 trong các bớc của 1
{
cho j=0 đến 15 trong các bớc của incr
{
i = j + ip
T = X(j)
X(i) = X(j) - T
X(j) = X(j) +T
nếu (iter không bằng 0) thì
{
cho k=1 đến ip-1 trong các bớc của 1
{
r = k*kk - 1
cho j=k đến 15 trong các bớc của 15
{
i=j+ip
T=X(i)*W(r)
X(i)=X(j)-1
X(j)=X(j)+T
}
}
}
kk=kk/2
ip= ip*2
inc=inc*2

int , int);
void main()
{
int i,k,m,N,n2,sign;
unsigned int *L;
float *wr,*wi,*xr,*xi;
char file_name[14];
FILE *fptr;
printf("\nEnter name of file containing data points-> ");
scanf("%s",file_name);
if((fptr=fopen(file_name,"rb"))==NULL)
{
printf("file %s does not exist.");
exit(1);
}
printf("Enter # of data points to be read >");
scanf("%d",&N);
m=(int)(log10((double)N)/log10((double)2.0));
k=1;
for(i=0;i<m;i++)
k<<=1 ;
if (k!=N)
{
printf("n Length of file has to be multiples of 2. ");
exit(1);
}
/* Allocating memory for bit reversal LUT. */
L=(unsigned int *)malloc(N*sizeof(unsigned int));
/* Generate Look-up table for bit reversal. */
bit_reversal(L,m,N);

Note: N=(2 to the power of m).
LUT will reside in LI]*/
{
unsigned int MASK,C,A,j,k,i;
for(k=0;k<N;k++)
{
MASK=1;
C=0;
for(i=0,j=m-1;i<m;i++,j )
{
A=(k&MASK)>>i;
A<<=j ;
C|=A;
MASK=MASK<<1;
}
L[k]=C;
}
}
void WTS(float *wr, float *wi, int N, int sign)
/* Generating LUT for twiddle factors.
Note:
sign=-1 for FFT, and
sign=1 for IFFT */
{
int n2,i ;
float theta;
91
n2=(N>>1)-1;
/* Generate look-up tables for twiddle factor. */
theta=2.0*pi/((float)N);

Ti=xi[i];
xr[i]=xr[j]-Tr;
xi[i]=xi[j]-Ti;
xr[j]=xr[j]+Tr;
xi[j]=xi[j]+Ti;
}
if(iter!=0)
{
for(k=1; k<ip; k++)
{
l=k*kk-1 ;
for(j=k; j<N; j+=incr)
92
{
i=j+ip ;
Tr=xr[i]*wr[l]-xi[i]*wi[l];
Ti=xr[i]*wi[l]+xi[i]*wr[l];
xr[i]=xr[j]-Tr;
xi[i]=xi[j]-Ti;
xr[j]=xr[j]+Tr;
xi[j]=xi[j]+Ti;
}
}
}
kk>>=1;
ip<<=1;
incr<<=1;
}
}
Chú ý rằng trong chơng trình 6.1 chúng ta giả thiết là dữ liệu đợc lu nh dãy của các ký

N
W
N
kfWkf
WkfWkfnF


=


=


=









++=
+=
12/
0
2/
1
2/

9!)]
2
()([)12(
)12(
2/
2/).12(
12/
0
+
+

=
++=+

nk
N
Nn
N
k
W
N
kfWkfnF

Chú ý rằng
0.1
2
==



=
++=

[ ]

])
2
()([)12(
2/
12/
0
kn
N
k
N
N
k
WW
N
kfkfnF


=
+=+

Đặt
f k f k f k
N

(6.39)

.)()12(
2/
12/
0
11
kn
N
N
k
WkfnF


=

=+
(6.40)
Các biểu thức (6.39) và (6.40) có thể biểu diễn bằng dới dạng biểu đồ bớm nh trong
hình 6.6.
Chúng ta có thể tiếp tục chia nhỏ các tổng cho trong các biểu thức (6.39) và (6.40),
tiếp tục làm nh vậy cho tới khi mỗi tổng giảm xuống chỉ còn lại một phần tử. Giải thuật
này giống nh giải thuật thuật toán phân chia thời gian và để lại cho bạn nh một bài tập
cho bạn. Một lu đồ cho FFT phân chia tần số với N = 4 trình bày trong hình 6.7. Bạn cần
chú ý đến bậc của dữ liệu đầu ra là bit đợc đảo. Phần mềm thực hiện thuật toán trên thì
rất giống phần mềm thực hiện FFT phân chia miền thời gian, và một chơng trình C đợc
cung cấp ở Chơng trình 6.2.
94
Có lẽ bạn sẽ tự hỏi: nếu phân chia miền thời gian đã thực hiện đợc công việc thì tại
sao lại phải xem xét thêm FFT phân chia tần số. Để trả lời câu hỏi này, chúng ta sẽ cần

incr=N;
for(iter=0; iter<m; iter++)
{
for(j=0; j<N; j+=incr)
{
i=i+ip;
Tr=xr[i];
Ti=xi[i];
xr[i]=xr[j]-Tr;
xi[i]=xi[j]-Ti;
xr[j]=xr[j]+Tr;
xi[j]=xi[j]+Ti;
}
for(k=1;k<ip; k++)
{
l=k*kk-1;
for(j=k; j<N; j+=incr)
96
0
1
2
3
4
5
6
7
0
2
4
6

5
13
3
11
7
15
0
8
4
12
2
10
6
14
1
9
5
13
3
11
7
15
0
4
8
12
2
6
12
14

6
7
8
9
10
11
12
13
14
15
n
W
8
−
n
W
4
−
n
W
2
−
n
W
−
{
i=j+ip;
Tr=xr[j]+xr[i];
Ti=xi[j]+xi[i];
diffr=xr[j]-xr[i];

ớm. Sơ đồ này giúp chúng ta thay đổi chơng trình 6.2 thành chơng trình 6.3.
Chơng trình 6.3 "FFTIP.C". Giảm lợc đầu vào FFT.
/****************************
* Program developed by: *
* M.A.Sid-Ahmed. *
* ver. 1.0 1992. *
* @ 1994 *
***************************/
/* FFT - input pruning routine. */
void bit_reversal(unsigned int *, int , int);
void WTS(float *, float *, int, int);
void FFTP(float *xr, float *xi, float *, float *, int,int,int, int);
void FFTP(float *xr, float *xi, float *wr, float *wi,
int m_output, int N_output,
int m_input, int N_input )
{
/* FFT pruning algorithm.
Deimation-in-frequency algorithm.
97
Note:
1. Noutput=2 to the power of m_output.
N_output=Number of samples in the output sequence.
M_input=Number of samples in the input sequence. This should also be a multiple of
2.
2. The output arrays are left in bit-reverse order.
You will need to use routine "bit-reversal" to place them in normal ascending order.
3. The twiddle factors are assumed to be stored in LUT's wr[I and wi[I. You will
need to use routine LUT for calculating and storing twiddle factors. */
int ip,k,kk,l,inc r,iter,j,i;
float Tr,Ti,diffr,diffi;

1!15
04812
26101415
91!!7111
5
08412210
6141951!
!11715
08412
21061419
51!!1171
5
W
-n
W
-6n
W
-4n
W
-2n
Hình 6.8 Lu đồ thuật toán giảm lợc đầu vào, N=4
for(iter=(m_output-m_input);iter<m_output;iter++)
{
for(j=0; j<N_output; j+=incr)
{
i=j+ip ;
Tr=xr [i];
Ti =xi [i] ;
xr[i]=xr[j]-Tr;
xi[i]=xi[j]-Ti;

Tính 1024 điểm trong phổ tần số dùng chơng trình giảm lợc đầu vào FFT.
99