96

l=k*kk-1 ;
for(j=k; j<N; j+=incr)
{
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 lưu như
dãy của các ký tự không dấu. Nếu bạn muốn xử lý trên một số dấu phẩy động bạn
cần thay đổi các câu lệnh mở và đọc dữ liệu trong file dữ liệu. Chương trình này
cũng cho phép lựa chọn FFT hoặc IFFT. Cho FFT, chương trình con "WTS( ) "
tính toán và lưu các hệ số dịch xoay trong một LUT được gọi lên vói tham số
"sign" được gán giá trị -1, ví dụ, WTS(wr,wi,N,-1) và cho IFFT,

N
W
N
kfWkf
WkfWkfnF




















12/
0
2/
1
2/

)]
2
()([)12(
)12(
2/
2/).12(
12/
0






nk
N
Nn
N
k
W
N
kfWkfnF

Chú ý rằng 0.1
2

 njnN

N
kfkfnF




 

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






(6.39)

98
.)()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à

Decimation-in-frequency algorithm.
Note :
1. N=2 to the power of m.

99

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[j and wiEj. You will need to use routine LUT
for calculating and storing twiddle factors. */

int ip,k,kk,l,incr,iter,j,i;
float Tr,Ti,diffr,diffi;

Hình 6.7 N = 4, phân chia miền tần số FFT.
ip= (N>>1) ;
kk=1;
incr=N;

for(iter=0; iter<m; iter++)
0
1
2
3
4
5
6
7

6
14
1
9
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
3
4
5
6
7
8
9
10
11
12
13
14
15
n
W
8n
W
4n
W
2n
W


100
{
for(j=0; j<N; j+=incr)
{
i=i+ip;
Tr=xr[i];

6.3.3 FFT giảm lược
Vấn đề có thể bắt đầu từ: cho 2
M
điểm dữ liệu, chúng ta sẽ phải làm thế nào
tính toán nhanh nhất khi dùng FFT có 2
L
điểm ra với M

L? Nếu M < L, có một
số bướm sẽ bị giảm lược (xem hình 6.8). Một giải thuật dựa trên tính toán các
phần tử của bướm khi việc tính toán tất cả các phép tính là không cần thiết trong

101
trường hợp M < L gọi là giải thuật giảm lược đầu vào FFT. Trong trường hợp M
> L thuật toán gọi là thuật toán giảm lược đầu ra FFT.

Thuật toán giảm lược đầu vào FFT. Trường hợp này sẽ làm hoàn thiện hơn
thuật toán phân chia tần số. Hình 6.8 giới thiệu trường hợp M = 1 và L = 4. Từ
hình 6.8 chúng ta nhận thấy (L-M) bước đầu tiên có các phần tử bướm và L bước
cuối cùng có toàn bộ các bướ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 *
***************************/

float Tr,Ti,diffr,diffi;
ip=N_output>>1;
kk=l ;
incr=N_output;
for(iter=0; iter<(m_output-m_input); iter++)
{
for(j=0; j<N_output; j+=incr)
{
i=i+ip ;
xr[i]=xr[j];
xi[i]=xi[j];
}
for(k=l; k<N_input; k++)
{
l =k*kk- 1 ;
for(j=k; j<N_output; j+=incr)
{
i=j+ip ;
xr[i]=xr[j]*wr[l]-xi[j]*wi[l];
xi[i]=xr[j]*wi[l]+xi[j]*wr[l];
}
}
kk<<=1;
ip>>=1;
incr>>=1;
}


11

13

15

0

4

8

12

2

6

10

14

1

5

9

13


13

3

11

7

15

0

8

4

12

2

10

6

14

1

9


Hình 6.8 Lưu đồ 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;

xr[j]=xr[j]+Tr;
xi[j]=xi[j]+Ti;
}
for(k=l; k<ip; k++)
{
l=k*kk-1;
for(j=k; j<N_output; j+=incr)
{
i=j+ip;

Tr=xr[j]+xr[i];
Ti=xi[j]+xi[i];

diffr=xr[j]-xr[i];
diffi=xi[j]-xi[i];
xr[j]=Tr;

xi[j]=Ti ;


Hình 6.9 Lưu đồ cho giảm lược đầu ra FFT, N = 4.

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

0
4

15

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

0
1

W

n=0 đến 7

W

n=0 đến 3


int m, int N, int m_output, int N_output)
{

/* FFT output pruning algorithm using
Decimation-in-time.
Note :
1. N=number of input samples
=2 to the power m.
N-output = number of output samples =2 to the power
motput.

2. The input arrays are assumed to be rearranged in
bit-reverse order.
You will need to use routine "bit-reversal" for that
purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[] and wi[]. You will need to use routine LUT
for calculating and storing twiddle factors.*/

int ip,k,kk,l,incr,iter,j,i;
float Tr,Ti;