Hình 6.4 Bước đầu tiên của lưu đồ FFT.
Hình 6.5 giới thiệu sơ đồ thuật toán FFT cho N = 16. Chú ý rằng do yêu cầu
ban đầu của chương trình mà dãy vào được sắp xếp lại và chứa ở X(k), ví dụ
X(k)
x(q)
k = 0 đến 15
Bạn sẽ chú ý trên sơ đồ rằng q là giá trị bit của k. Cho N = 24 = 16 chúng ta phải có bốn bước trong lưu đồ. Trong mỗi bước cần phải có tám bướm. Trong mỗi bướm chỉ có một phép nhân phức, hai phép cộng hoặc trừ phức. Tổng số phép nhân phức là 8/2 . 4. Tổng quát cho N = 2r số phép nhân phức là (N/2) . r = (N/2 ) log2 N và số phép cộng là Nlog2N. 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 W0 = 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 = 210. 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ừ lưu đồ. Tuy nhiên, chúng ta có thể nghiên cứu lưu đồ và vạch ra các bước có thể dùng để phát triển một chương trình. Từ lưu đồ 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)
W0 = 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;
}
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)
86
{ 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 W4 = W(3). Chú ý rằng chúng ta coi rằng các hệ
số xoay W, W2 , ... , W7 đã được tính và được chứa trong W(0), W(1), ...W(6).
for (j=0 đến 15 với bước tăng 4) { T=X(j)W(3);X(j+2)=X(j) - T; X(j)=X(j) + T; }
Bước thứ ba. Chúng ta có : 1. Hai bướm với trọng lượng bằng 1.
for (j=0 đến 15 với bước tăng 8) { T=X(j); X(j+4)=X(j) - T; X(j)=X(j) + T; }
87
nW 8
nW 4
nW 2
nW
n=0 n=0 ®Õn 1 n=0 ®Õn 3 n=0 ®Õn 7
0 0 0 0 0000 0 0000
8 4 2 1 0001 8 1000 0 0 4 8 4 2 0010 4 0100
4 12 12 6 3 0011 12 1100 0 2 2 8 4 0100 2 0010 0
10 6 10 5 0101 10 1010 2 0 6 12 12 6 0110 6 0110 4 0 14 14 14 7 0111 14 1110 6 0 4 1 1 1 8 1000 1 0001 0
9 5 3 9 1001 9 1001 1 0 5 9 5 10 1010 5 0101 2 0 13 13 7 11 1011 13 1101 3 0 4 3 3 9 12 1100 3 0011 4 0 11 7 11 13 1101 11 1011 5 2 0 7 11 13 14 1110 7 0111 6 4 0 15 15 15 15 1111 15 1111 7 6 0 4
bíc 0 bíc 1 bíc 2 bíc 3
BËc cña d·y vµo biÓu diÔn d¹ng nhÞ ph©n BËc cña d·y ra biÓu diÔn d¹ng nhÞ ph©n
Hình 6.5 Lưu đồ thuật toán thuật toán phân chia miền thời gian.
2. Hai bướm với trọng lượng bằng W (1) = W2.
for (j=1 đến 15 với bước tăng 8) { T=X(j)W(1); X(j+4)=X(j) - T; X(j)=X(j) + T; }
88
3. Hai bướm với trọng lượng bằng W(3) = W4.
for (j=2 đến 15 với bước tăng 8) { T=X(j)W(3); X(j+4)=X(j) - T; X(j)=X(j) + T; }
f
f
F(2n)
F(2n)
)(10 k
)(10 k
F(2n+1)
F(2n+1)
f
f
)(11 k
)(11 k
(a)
(b)
Hình 6.6 (a) Bướm cho thuật toán phân chia miền tần số;(b) Một bướm đơn giản.
4. Hai bướm với trọng lượng bằng W (5) = W6 .
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
89
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) = W2. T = X(1)W(1)
X(2+8)= X(0) - T X(2) = X(2) +T . . .
8. Một bướm với trọng lượng bằng W (6) = W7.
T = X(7)W(6) X(7+8)= X(7)-T X(7) = X(7) +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 {
90
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 } Thuật toán trên có thể dễ dàng mở rộng cho tất cả các trường hợp của N. Chỉ có một lĩnh vực còn lại cần phải giải thích là sự sắp xếp lại các dãy dữ liệu đầu vào. Điều này có thể tạo ra dễ dàng nếu chúng ta tạo ra một bảng (LUT) L(i), L(i) là các giá trị đảo ngược bit của i. Nếu dữ liệu được đọc từ một file thì tất cả các việc mà chúng ta phải làm là chuyển địa chỉ vùng của chúng trong file qua bảng LUT và lưu các dữ liệu này trong địa chỉ chứa kết quả trong dãy đầu vào, X. Bước này có thể chuyển sang ngôn ngữ C như sau:
for (i=0; i
Chương trình 6.1 “FFTDT.C” FFT 1-D Thập phân trong miền thời gian.
/***********************************
* Program developed by: *
* M.A.Sid-Ahmed. *
* ver. 1.0 1992. *
* @ 1994 *
*********************************/
/* FFT - Decimation-in-time routine with examplemain
programing proper usage. */
#define pi 3.141592654
91
int , int);
#include
#include
#include
#include
void bit_reversal(unsigned int *, int , int);
void WTS(float *, float *, int, int);
void FFT(float *xr, float *xi, float *, float *,
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
92
/* Generate Look-up table for bit reversal. */
bit_reversal(L,m,N);
/* Allocating memory for FFT arrays ( real and imag.)
*/
xr=(float *)malloc(N*sizeof(float));
xi=(float *)malloc(N*sizeof(float));
/* Setting-up the data in real and imag. arrays.*/
for(i=0;i>1)-1;
wr=(float *)malloc(n2*sizeof(float));
wi=(float *)malloc(n2*sizeof(float));
/*Generating LUT for
twiddle factors. */
WTS(wr,wi,N,-1);
/* Taking FFT. */
FFT(xr, xi, wr, wi, m, N);
printf("Enter file name for storing FFT output.--->");
scanf("%s",file_name);
fptr=fopen(file_name,"w");
for(i=0;i
93
unsigned int MASK,C,A,j,k,i;
for(k=0;k>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;
n2=(N>>1)-1;
/* Generate look-up tables for twiddle factor. */
theta=2.0*pi/((float)N);
for(i=0;i
94
You will need to use routine "bitreversal" for
{
/* FFT algorithm,
Decimation-in-time algorithm.
Note:
1. N=2 to the power of m.
2. The input arrays are assumed to be rearranged in
bit-reverse order.
that purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[I and wi[j. 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;
ip=1;
kk=(N>>1);
incr=2 ;
for(iter=0; iter
95
Chương trình 6.1 “FFTDT.C” FFT 1-D Thập phân trong miền thời gian.
/*********************************** * Program developed by: * * M.A.Sid-Ahmed. * * ver. 1.0 1992. * * @ 1994 * *********************************/ /* FFT - Decimation-in-time routine with examplemain programing proper usage. */ #define pi 3.141592654
91
int , int);
#include
#include
#include
#include
void bit_reversal(unsigned int *, int , int);
void WTS(float *, float *, int, int);
void FFT(float *xr, float *xi, float *, float *,
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
92
/* Generate Look-up table for bit reversal. */
bit_reversal(L,m,N);
/* Allocating memory for FFT arrays ( real and imag.)
*/
xr=(float *)malloc(N*sizeof(float));
xi=(float *)malloc(N*sizeof(float));
/* Setting-up the data in real and imag. arrays.*/
for(i=0;i>1)-1;
wr=(float *)malloc(n2*sizeof(float));
wi=(float *)malloc(n2*sizeof(float));
/*Generating LUT for
twiddle factors. */
WTS(wr,wi,N,-1);
/* Taking FFT. */
FFT(xr, xi, wr, wi, m, N);
printf("Enter file name for storing FFT output.--->");
scanf("%s",file_name);
fptr=fopen(file_name,"w");
for(i=0;i
93
unsigned int MASK,C,A,j,k,i;
for(k=0;k>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;
n2=(N>>1)-1;
/* Generate look-up tables for twiddle factor. */
theta=2.0*pi/((float)N);
for(i=0;i
94
You will need to use routine "bitreversal" for
{
/* FFT algorithm,
Decimation-in-time algorithm.
Note:
1. N=2 to the power of m.
2. The input arrays are assumed to be rearranged in
bit-reverse order.
that purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[I and wi[j. 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;
ip=1;
kk=(N>>1);
incr=2 ;
for(iter=0; iter
95
92
/* Generate Look-up table for bit reversal. */
bit_reversal(L,m,N);
/* Allocating memory for FFT arrays ( real and imag.)
*/
xr=(float *)malloc(N*sizeof(float));
xi=(float *)malloc(N*sizeof(float));
/* Setting-up the data in real and imag. arrays.*/
for(i=0;i>1)-1;
wr=(float *)malloc(n2*sizeof(float));
wi=(float *)malloc(n2*sizeof(float));
/*Generating LUT for
twiddle factors. */
WTS(wr,wi,N,-1);
/* Taking FFT. */
FFT(xr, xi, wr, wi, m, N);
printf("Enter file name for storing FFT output.--->");
scanf("%s",file_name);
fptr=fopen(file_name,"w");
for(i=0;i
93
unsigned int MASK,C,A,j,k,i;
for(k=0;k>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;
n2=(N>>1)-1;
/* Generate look-up tables for twiddle factor. */
theta=2.0*pi/((float)N);
for(i=0;i
94
You will need to use routine "bitreversal" for
{
/* FFT algorithm,
Decimation-in-time algorithm.
Note:
1. N=2 to the power of m.
2. The input arrays are assumed to be rearranged in
bit-reverse order.
that purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[I and wi[j. 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;
ip=1;
kk=(N>>1);
incr=2 ;
for(iter=0; iter
95
93
unsigned int MASK,C,A,j,k,i;
for(k=0;k>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;
n2=(N>>1)-1;
/* Generate look-up tables for twiddle factor. */
theta=2.0*pi/((float)N);
for(i=0;i
94
You will need to use routine "bitreversal" for
{
/* FFT algorithm,
Decimation-in-time algorithm.
Note:
1. N=2 to the power of m.
2. The input arrays are assumed to be rearranged in
bit-reverse order.
that purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[I and wi[j. 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;
ip=1;
kk=(N>>1);
incr=2 ;
for(iter=0; iter
95
94
You will need to use routine "bitreversal" for
{
/* FFT algorithm,
Decimation-in-time algorithm.
Note:
1. N=2 to the power of m.
2. The input arrays are assumed to be rearranged in
bit-reverse order.
that purpose.
3. The twiddle factors are assumed to be stored in
LUT's wr[I and wi[j. 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;
ip=1;
kk=(N>>1);
incr=2 ;
for(iter=0; iter
95
95
