Giáo trình Xử lý ảnh - Chương 6: Biến đổi Fourier rời rạc

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: (6.1) Nếu h(k1,k2) chỉ có chỉ tồn tại với k1 ³ 0, k 2 ³ 0 và tổng quát được xác định trong miền hữu hạn có kích thước N ´ N thì (6.2) Công thức này chứng tỏ rằnglà tuần hoàn, chu kỳ tuần hoàn là 2p. Nếu chúng ta lấy mẫu dưới dạng w1, w2, và miền xác định là (0 Ê w1 Ê 2p) và (0 Ê w2 Ê 2p), N ´ N mẫu, chúng ta có thể viết: và (6.3) vì thế (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 : (6.5) Công thức này có thể viết lại dưới dạng (6.6) ở đây f(k) = f(kT) và WN = e- j2 /N . WN được gọi là hạt nhân của phép biến đổi. Tổng quát, F(n) có dạng (6.7) Ký hiệu A(n), f (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 (6.8) Chứng minh: Từ định nghĩa của DFT (6.9) Đặt Nếu (k = m) thì S = N. Nếu (k ạ m), chúng ta có thể viết: S = 1 + WN (k -m ) + WN 2(k -m ) + ... + WN (N-1)(k -m ) hoặc Khi e j2p(k-m) = 1 và e j2p/N. (k-m) ạ 1 với (k ạ m), vì vậy S = 0 với (k ạ m ). Vì vậy, biểu thức (6. 9) có thể rút gọn thành 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 Chú ý từ biểu thức (6.8) ta có thể dễ dàng chứng minh: (6.10) Mặc dù f(k) được xác định trên miền k ẻ [0,N], nó vẫn là tín hiệu tuần hoàn với chu kỳ NT. (T được bao hàm và rút ra từ biểu thức 6.5). 6.2.2 Một vài tính chất của DFT Tuyến tính. Nếu ta có hai dãy tuần hoàn cùng f1(n) và f2(n), và cả hai dãy này tuần hoàn với chu kỳ N, được dùng để tính f3(k) = af1(k) + bf2(k) (6.11) là kết quả của biến đổi DFT f3(n) cho bởi F3(n) = aF1(n) + bF2(n) (6.12) ở đây a, b là hằng số và F1(n) = DFT của f1(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. Nếu f(k) là thực thì (6.13) Dấu * có nghĩa là liên hợp phức. Tích chập tuần hoàn. Coi f1(k) và f2(k) là hai dãy tuần hoàn có chu kỳ N, với biến đổi Fourier rời rạc là F1(n) và F2(n). Xem xét tích F(n1).F(n2) khi và tại các vị trí n1 = n2 = n Đặt f3(k) = IDFT của F1(n).F2(n) hay vì vậy cho k = k1 + k2 + lN các trường hợp còn lại. Chú ý là ở đây l là số nguyên. Vì vậy mà (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. Bài tập 6.1 Cho hai dãy sau 0 Ê k1 Ê4 các trường hợp còn lại. 0 Ê k1 Ê 5 các trường hợp còn lại. Tính bằng tay tích chập của hai dãy trên. Vẽ một lưu đồ biểu diễn thuật toán. Làm lại phần 1, nhưng lần này sử dụng tích chập tuần hoàn. 3. Lập một chương trình C rút ra f3(k) từ biểu thức f3(k) = IDFT{DFT[f1(k)]. DFT[f1(k)]}. So sánh kết quả của phần 1 và phần hai. 4. Bây giờ thêm các số không vào f1(k) và f2(k) để chu kỳ của chúng = 5 + 6 - 1. Làm lại phần 3 và so sánh kết quả. 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ì N2 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 = 210 ta sẽ phải tính 220 = 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= 2r (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, Chia dãy f(k) thành hai dãy, một dãy được rút ra từ phần tử chẵn và một dãy từ những phần tử lẻ. Đó là, k chẵn k lẻ Chúng có thể viết thành (6.15) Chú ý là vì thế đặt Ta được (6. 16) Đặt (6.17) (6.18) Viết lại biểu thức (6.16) chúng ta được (6.19) Cũng như vậy, phát triển cho một biểu thức (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 W16-n thường gọi là trọng lượng hay hệ số xoay. Hai biểu thức này biểu diễn bước cuối cùng trong lưu đồ tính toán của hình 6.2. Bây giờ xem xét biểu thức Xử lý như trên chúng ta có Dễ thấy đặt F10(n) F(n) F11(n) F(n+8) F10(n) F(n) F11(n) F(n+8) 1 W16-n n Hình 6.1 (a) Bướm; (b) Biểu diễn rút gọn. 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 F10(n) F(n) F11(n) X(k) X(k) Hình 6.2 Bước cuối cùng của thuật toán biến đổi FFT phân chia miền thời gian. X(k) ký hiệu vector chứa giá trị được tính qua phép biến đổi FFT. Vì vậy, (6.21) (6.22) ở đây (6.23) (6.24) Tương tự (6.25) (6.26) ở đây (6.27) (6.28) và Biểu thức (6.21), (6.22), (6.25) và (6.26) có thể biểu diễn bằng sơ đồ hình 6.3. Biểu thức (6.23), (6.24), (6.27) và (6.28) có thể tiếp tục chia nhỏ ra như các bước đã làm ở trên như sau: (6.29) (6.30) (6.31) (6.32) (6.33) (6.34) (6.35) (6.36) ở đây (6.37) (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 lưu đồ hình 6.4.Mỗi phần tử từ F30(n) đến F37(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 lưu đồ. 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 0 F20(n) 2 F21(n) 4 6 F22(n) 0 F23(n) 4 2 6 F10(n) F11(n) X(k) X(k) W-2n Hệ số xoay n=0 đến 3 Hình 6.3 Bước thứ hai sau bước cuối cùng trong thuật toán FFT. 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 F30(n) F31(n) F32(n) F33(n) F34(n) F35(n) F36(n) F37(n) X(k) X(k) Dãy đầu vào đã được sắp xếp lại 0 0 0 0 0 0 0 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ụ 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) { 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; } 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; } 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; } 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 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 { 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<N; i++) fscanf (fptr, “ %f ”, &X[L[i]]); Kết quả của LUT được chuyển thẳng và được cung cấp với chương trình của thuật toán tính FFT trong Listing 6.1 dưới dạng modun con dưới tên “ bit_reversal(...)”. 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 #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 *, 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); /* 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<N;i++) { k=L[i] ; xr[k]=(float)getc(fptr); xi[k]=0.0 ; } fclose(fptr); /* Allocating memory for twiddle factors. */ n2=(N>>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<N;i++) fprintf(fptr," %e %e",xr[i], xi[i]); fclose(fptr); } void bit_reversal(unsigned int *L, int m, int N) /* Routine for generating LUT for bit reversal. 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; n2=(N>>1)-1; /* Generate look-up tables for twiddle factor. */ theta=2.0*pi/((float)N); for(i=0;i<n2;i++) { wr[i]=(float)cos((double)((i+1)*theta)); wi[i]=(float)sin((double)((i+1)*theta)); if(sign==(int)(-1)); wi[i]=-wi[i]; } } /********************************************/ void FFT (float *xr, float *xi, float *wr, float *wi, int m, int N) { /* 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. You will need to use routine "bitreversal" for 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<m; iter++) { for(j=0; j<N; 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; } if(iter!=0) { for(k=1; k<ip; k++) { 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, WTS(wr,wi,N,1). Với IFFT, bạn cần chia dãy ra cho N trong chương trình gọi hoặc là chương trình chính. Bài tập 6.2 Kiểm tra chương trình FFT bằng cách làm lại chương trình 6.1. Chú ý rằng trong trường hợp này bạn phải thêm các giá trị 0 để làm cho các dãy có chiều dài 24 = 16 và tất nhiên là lớn hơn chiều dài dãy nhỏ nhất đòi hỏi là (6 + 5 - 1). Mối tương quan của hai dãy cho kết quả trong một tín hiệu tuần hoàn có chu kỳ bằng 16. 6.3.2 Thuật toán phân chia tần số. Thay vì chia dãy vào thành các vị trí chẵn và lẻ, chúng ta sẽ đưa ra một chương trình giống như chương trình trên nhưng lần này ta bắt đầu từ dãy ra. Chương trình này bao gồm các bước sau: Bây giờ, chia dãy F(n) thành hai dãy dựa trên giá trị chẵn và lẻ của n. Chú ý rằng Vì vậy Đặt Vì vậy (6.39) (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 lưu đồ 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. 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 xem xét phần kế tiếp, FFT giảm lược. Chương trình 6.2 “FFTDF” FFT phân chia tần số. /**************************** * Program developed by: * * M.A.Sid-Ahmed. * * ver. 1.0 1992.1994 * *****************************/ /* FFT - Decimation-in-frequency routine.*/ #define pi 3.141592654 void bit_reversal (unsigned int *, int, int); void WTS(float *, float *, int, int) ; void FFT(float *xr, float *xi , float, float, int, int); void FFT (float *xr, float *xi, float *wr, float *wi, int m, int N) { /* FFT algorithm. Decimation-in-frequency algorithm. Note : 1. N=2 to the power of m. 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++) { 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) { 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 ; Tr=diffr*wr[1]-diffi*wi[1]; Ti=diffr*wi[1]+diffi*wr[1]; xr[i]=Tr; xi[i]=Ti ; } } kk<<=1; ip>>=1 ; incr>>=1; } } 6.3.3 FFT giảm lược Vấn đề có thể bắt đầu từ: cho 2M đ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ó 2L điểm ra với M ạ L? Nếu 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 * ***************************/ /* 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. 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; 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; } 0 1 0 2 4 6 8 10 12 14 1 3 5 7 9 11 13 15 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15 0 8 4 12 2 10 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 W-n W-6n W-4n W-2n 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 ; Tr=diffr*wr[l]-diffi*wi[l]; Ti=diffr*wi[l]+diffi*wr[l]; xr[i]=Tr; xi[i]=Ti; } } kk<<=l; ip>>=l; incr>>=l; } } Bài tập 6.3 1.Cho dãy đầu vào : x(k) = 1 k = 0,1,2, ... , 31. x(k) = 0 các trường hợp còn lại. Tính 1024 điểm trong phổ tần số dùng chương trình giảm lược đầu vào FFT. 2. Thêm các giá trị 0 vào dãy để làm cho chiều dài dãy thành 1024. Bây giờ tính FFT scủa dãy dùng chương trình FFT phân chia tần số không giảm lược. So sánh thời gian xử lý của phần 1 và 2. Thuật toán FFT giảm lược đầu ra. Giải thuật phân chia miền thời gian thì thích hợp cho thuật toán giảm lược đầu ra hơn là giải thuật phân chia miền tần số. Lý do là đầu ra trong giải thuật phân chia miền thời gian không phải sắp xếp lại. Hình 6.9 giới thiệu trường hợp với M=4 và L=1. 0 2 4 6 8 10 12 14 1 3 5 7 9 11 13 15 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15 0 8 4 12 2 10 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 1 W n=0 đến 7 W n=0 đến 3 W n=0 đến 1 W n=0 Hình 6.9 Lưu đồ cho giảm lược đầu ra FFT, N = 4. Chương trình 6.4 “FFTOP.C” Giảm lược đầu ra FFT. /***************************** * Program developed by: * * M.A.Sid-Ahmed. * * ver. 1.0 1992. * * @ 1994 * *****************************/ /* FFT - output pruning using Decimation-in-time routine. */ # define pi 3.141592654 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, 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; ip=1; kk=(N>>1 ); incr=2; for(iter=0; iter<m_output; iter++) { for(j=0; j<N; j+=incr) { i=j+ip; Tr=xr[i]; Ti=xi[i]; xr[i]=xr[i]-Tr; xi[i]=xi[i]-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) { 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; } for(iter=m_output; iter<m; iter++) { for(j=0; j<N; j+=incr) { i=j+ip; xr[j]=xr[j]+xr[i]; xi[j]=xi[j]+xi[i]; } for(k=1; k<N_output; k++) { 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[i]; xr[j]=xr[j]+Tr; xi[j]=xi[j]+Ti; } } kk>>=1; ip<<=1; incr<<=1; } } Bài tập 6.4 Cho một bộ lọc thông cao Butterword: ở đây w0 = 0.3p, phát triển một chương trình C tính ra 1024 mẫu của đặc tuyến tần số trên các khoảng bằng nhau của tần số w từ 0 đến p. Dùng các mẫu này tính năm mẫu đầu tiên của đáp ứng xung dùng FFT giảm lược. 6.4 2-D FFT Một DFT hai chiều của tín hiệu lấy mẫu hai chiều h(k1,k2) cho bởi (6.41) ở đây n1 = 0,1,2 , ..., N-1 n2 = 0,1,2,..., N-1 Biểu thức trong hai dấu tổng gọi là hạt nhân của phép biến đổi. H(n1,n2), trong trường hợp tổng quát, đầy đủ có thể biểu diễn theo: Trong không gian ba chiều, A(n1,n2) và f(n1,n2) nằm tại vị trí của n1 và n2 và gọi là phổ tần số và phổ pha của H(n1,n2). 6.4.1 Biến đổi ngược 2-D DFT Hàm h(k1,k2) là biến đổi ngược của 2-D DFT (IFFT) của H(n1,n2) và được cho bởi biểu thức (6.42) 6.4.2 Một số tính chất của 2-D DFT Chuyển đổi. Từ định nghĩa của 2-D DFT và IDFT cho thấy (6.43) (6.44) Điều đó có nghĩa là một dịch chuyển pha tuyến tính trong một miền biểu diễn bằng một dịch chuyển hằng số trong một miền khác. Xem xét biểu thức (6.43), trường hợp đặc biệt khi a = b = N/2. Hay là (6.45) Nói cách khác, bằng cách nhân vào mỗi điểm trước khi lấy DFT, chúng ta sẽ rút ra được một phổ tần số mà điểm tần số (0,0) của nó sẽ nằm giữa mảng 2-D. Biểu thức này rất hữu dụng trong hiển thị phổ tần số, phổ biên độ và lọc dùng DFT. Từ biểu thức (6.44) chúng ta rút ra kết luận rằng dịch chuyển một hằng số trong ảnh sẽ không tác động đến phổ biên độ. (6.46) Biểu thức (6.46) cũng quan hệ đến bộ lọc 2-D. Xem xét đặc tính của bộ lọc 2-D cho bởi ở đây A(n1,n2) là phổ biên độ. Nếu một ảnh với phổ tần số cho bởi I(n1,n2) được lọc qua bộ lọc có đặc tuyến pha tuyến tính cho bởi biểu thức ở trên, kết quả sẽ là (6.47) ở đây if (n1-a, n2-b) ký hiệu cho ảnh đã được lọc. Một bộ lọc với đặc tuyến pha tuyến tính có nghĩa là không dịch chuyển biên độ. Trong khi đó nếu bộ lọc có đặc tuyến pha không tuyến tính thì pha của ảnh cũng bị biến dạng. Lý do của sự biến dạng này là tất cả các điểm đều phải chịu một sự dịch chuyển vị trí khác nhau tuỳ theo vị trí của ảnh. Tổng quát, ảnh đã được lọc có thể cho bởi ở đây f là hàm dịch chuyển vị trí. Chú ý rằng một ảnh biến dạng pha sẽ xuất hiện trên màn hình như một ảnh mờ . Tính đối xứng liên hợp và tuần hoàn. Biến đổi2-D DFT và IDFT tuần hoàn với chu kỳ N có nghĩa là : (6.48) và (6.49) Biến đổi DFT đối xứng liên hợp khi (6.50) hoặc (6.51) Quay. Nếu chúng ta đặt k1 và k2 dưới dạng thì Và tương tự cho n1, n2 hoặc từ định nghĩa của DFT chúng ta có thể có (6.52) Điều đó có nghĩa là nếu ảnh bị quay đi một góc f0 thì phổ của nó cũng bị quay đi một góc như vậy. Phân phối và chia độ. Từ biểu thức (6.1) chúng ta dễ thấy (6.53) ở đây a,b là những độ chia. Cũng như vậy (6.54) 6.4.3 Giá trị trung bình Mức cường độ sáng trung bình trong một ảnh cho bởi : (6.55) hoặc từ biểu thức (6.1) ta có thể viết Điều này có nghĩa là H(0,0) biểu diễn mức sáng của ảnh. 6.4.4 Tích chập và sự tương quan Tích chập của tín hiệu 2-D h1(k1,k2) và h2(k1,k2) cho bởi (6.56) Nếu h1(k1,k2) được xác định trên miền và nếu h2(k1,k2) được xác định trên miền thì chúng ta có thể thấy rằng nếu hai tín hiệu có giá trị zero ngoài miền xác định của chúng thì M = A + C - 1 và N = B + D - 1. Tích chập của hai tín hiệu 2-D được viết trong dạng ký hiệu như sau: Có thể thấy rằng (6.57) Điều này có nghĩa là, tích chập trong miền không gian biến thành phép nhân bình thường trong miền tần số. Tính chất này có thể dùng cho lọc 2-D qua DFT. Chúng ta cần nhớ lại rằng bạn đã dùng kỹ thuật lọc FIR trong các chương trước cho chức năng này. Khi áp dụng các lọc bộ lọc FIR cho chức năng lọc bạn cần lấy tín hiệu khoảng cách 2-D đã được biến thành tín hiệu có chu kỳ trước khi tiến hành lấy DFT. Sự không đồng bộ của chu kỳ trong biến đổi này cũng gây ra lỗi như trong biến đổi 1-D. Vì vậy, để tránh trường hợp này ta cần thêm các số 0 vào cả hai các hàm không gian để cho chúng có kích thước M ´ N với M ³ A + C - 1 và N ³ B + D - 1. Tương quan hoặc tương quan chéo của tín hiệu 2-D định nghĩa bởi (6.58) Biểu thức này được viết dưới dạng ký hiệu (6.59) Tương quan chéo thường được gọi là lọc kết hợp và dùng để phát hiện ra phần đầu dấu hiệu các vết sắc nổi trên ảnh. Nó có thể cho thấy rằng (6.60) 6.5 2-D FFT Biểu thức cho 2-D DFT có thể dưới dạng (6.61) Đặt (6.62) với k1 = 0,1,2,..., N-1. Vì vậy mà biểu thức (6.61) có thể viết lại dưới dạng . (6.63) với n1=0,1,2,...,N-1. Biểu thức (6.62) có thể mở rộng ra thành . . . ...v.v... Mỗi biểu thức trên biểu diễn DFT của một hàng trong ảnh. Biểu thức (6.63) cũng có thể mở rộng ra thành: . . . ...v.v... Các biểu thức này dẫn ta đến giải thuật sau đây tính FFT hai chiều: 1. Tính 1-D FFT cho tất cả các hàng, và chứa kết quả vào mảng trung gian. 2. Dịch chuyển mảng trung gian. 3. Rút ra 1-D FFT cho tất cả các hàng của mảng dịch chuyển trung gian. Kết quả là dịch chuyển của mảng 2-D FFT. Chúng ta có thể viết biểu thức (6.61) có dạng (6.64) Nếu chúng ta đặt (6.65) với k2 = 0,1,2,..., (N-1) thì (6.66) với n1 = 0,1,2,..., (N-1). Các biểu thức này dẫn chúng ta đến thuật toán tính 2-D FFT sau 1. Dịch chuyển file ảnh. 2. Tính FFT theo từng hàng một của ảnh đã được đọc. 3.Dịch chuyển kết quả trung gian. 4. Rút ra một hàng kề hàng FFT của dịch chuyển kết quả trung gian. Kết quả ta sẽ được 2-D FFT. Trong cả hai phương pháp dùng để rút ra 2-D FFT, kết quả trung gian đều đã phải dịch chuyển. Phương pháp đầu tiên thường hay được sử dụng hơn vì nó chỉ yêu cầu một phép toán dịch chuyển. Kết quả là một dịch chuyển của mảng 2-D FFT, có thể dùng trực tiếp dưới dạng ấy mà không đòi hỏi một phép dịch chuyển thứ hai. Chắc bạn sẽ có một câu hỏi rằng tại sao chúng ta cần phải dịch chuyển. Lý do của sự dịch chuyển này là hệ thống của bạn có thể không có đủ bộ nhớ kích hoạt (RAM) để lưu trữ kết quả trung gian hoặc là FFT của ảnh. Nếu bạn có đủ bộ nhớ RAM thì việc dịch chuyển này là không cần thiết, và bạn có thể đọc thẳng từng cột từ bộ nhớ kích hoạt. Dù sao đi chăng nữa thì sự lựa chọn vẫn là đọc thẳng từng cột và có kết quả trung gian chứa theo hàng. Nếu là như vậy, chúng ta sẽ cần N ´ N dữ liệu thêm vào từ đĩa cứng, yêu cầu thời gian nhiều hơn. Nói một cách khác, dịch chuyển file dẫn đến từng hàng một trong FFT của kết quả trung gian, đòi hỏi nhiều hơn N lần truy nhập đĩa. Câu hỏi bây giờ là làm thế nào chúng ta có thể dịch chuyển một file trong trường hợp không thể chuyển tất cả dữ liệu một lần vào bộ nhớ kích hoạt. Trong phần tiếp theo chúng ta sẽ đề cập đến phương pháp Eklundh để giải quyết vấn đề này. 6.5.1 Ma trận dịch chuyển từ bộ nhớ ngoài Thuật toán được giải thích rõ ràng nhất bằng một ví dụ đặc biệt. Xem xét ma trận có kích thước 7 ´ 7 ở hình 6.10. Các bước của thuật toán thể hiện rõ ràng trên hình 6.10. Bạn cần chú ý rằng chương trình đòi hỏi ba lần lặp lại để dịch chuyển một ma trận có kích thước 23 ´ 23. Trong tất cả các lần lặp lại bạn cần phải giữ lại trong bộ nhớ kích hoạt tại hai hàng cuối cùng, cho phép lặp, tất cả là yêu cầu yêu cầu N lần truy nhập đĩa cho xử lý một ảnh có kích thước N ´ N. Nếu N = 2r thì r ´ N số lần truy nhập đĩa để dịch chuyển một ảnh, ít hơn nhiều so với N ´ N lần truy nhập trong cách xử lý đọc một ảnh cơ bản từng khối một. Số lần truy nhập đĩa có thể giảm xuống bởi đọc, ví dụ, bốn hàng hoặc tám hàng một lúc. Thuật toán trong trường hợp tổng quát có thể coi như là sự phát triển của giải thuật FFT. Bạn cần phát triển hai lưu đồ, một để lựa chọn hàng có thể đã được thêm dữ liệu vào và một để phát hiện ra phần tử đã thay đổi. Những thuật toán này coi như là một bài tập. Mã của chương trình nguồn cho tất cả các thuật toán này cho ở chương trình 6.5. Chương trình 6.5 “TRANSPOS.C” chương trình cho dịch chuyển một ma trận. /****************************** * Program developed by: * * M.A.Sid-Ahmed. * * ver. 1.0 1992. * * @ 1994 * ******************************/ /*This program is for obtaining the transpose of a large binary file stored on secondary memory. We assume that the file is large to the point that it would not fit in active memory, and the data type in the file is character.*/ #include #include #include #include #include #include void transpose(FILE *, int, int); void main() { int N,n; char file_name[11]; FILE *fptr; float nsq; clrscr(); printf("Enter file-name to be transposed -->"); scanf("%s",file_name); fptr=fopen(file_name,"rb+"); if(fptr==NULL) { printf("%s does not exist. "); exit(1); } nsq=filelength(fileno(fptr)); N=sqrt((double)nsq); n=(int)(log10((double)(N))/log10((double)(2.0))); clrscr(); transpose(fptr,N,n); fclose(fptr); } void transpose(FILE *fptr, int N, int n) /* Algorithm */ { int N1,inc; int iter,i,k; int k1,inc1; int k2,j,k3,k4; char *buff1,*buff2,tmp; long loc; buff1=(char *)malloc(N*sizeof(char)); buff2=(char *)malloc(N*sizeof(char)); N1=N/2; inc=1; inc1=2; Hình 6.10 Thuật toán của Eklundh cho dịch chuyển một ma trận. for(iter=0;iter<n;iter++) { gotoxy(1,2); printf("iteration # %4d",iter+1); k1=0 ; for(k=0;k<N1;k++) { for(i=k1;i<(k1+inc);i++) loc=(long)(N)*(long)(i); if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1); } else { gotoxy(1,3); printf("Reading row # %4d",i); for(k4=0;k4<N;k4++) *(buff1+k4)=fgetc(fptr); { loc=(long)(N)*(long)(i+inc); } if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1) ; } else { gotoxy(1,4); printf("Reading row # %4d",i+inc); for(k4=0;k4<N;k4++) *(buff2+k4)=fgetc(fptr); } k3=0; for(k2=0;k2<N1;k2++) { for(j=k3;j<(k3+inc);j++) { tmp=*(buff1+j+inc); *(buff1+j+inc)=*(buff2+j); *(buff2+j)=tmp; } k3+=inc1 ; } loc=(long)(N)*(long)i; if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit( 1 ) ; } else { gotoxy(1,3); printf("writing row # %4d",i); for(k4=0;k4<N;k4++) putc((char)(*(buff1+k4)),fptr); } loc=(long)(N)*(long)(i+inc); if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1) ; } else { gotoxy(1,4); printf("writing row # %4d",i+inc); for(k4=0;k4<N;k4++) putc((char)(*(buff2+k4)),fptr); } } k1+=inc1 ; } inc*=2; inc1*=2; N1/=2; } Để kiểm tra chương trình 6.5 chúng ta sẽ áp dụng thuật toán này lên ảnh cho trong hình 6.11. ảnh này chứa trên đĩa đi kèm theo cuốn sách này dưới file có tên là “MOHSEN.IMG”. Chương trình 6.6 “ FFT2D.C” 2-D FFT /****************************** * Program developed by: * * M.A.Sid-Ahmed. * * ver. 1.0 1992. * * @ 1994 * ******************************/ /* 2D-FFT - Using Decimation-in-time routine.*/ #define pi 3.141592654 #include #include #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 *, int, int ) ; void transpose(FILE *, int, int); void FFT2D(FILE *, FILE *, float *, float*, unsigned int *,int,int,int); Hình 6.11 ảnh đã được dịch chuyển, "MOHSEN.IMG". void main() { int N,n2,m,k,i; unsigned int *L; float *wr , *wi; char file_name[14]; FILE *fptr,*fptro; double nsq; clrscr(); printf(" Enter name of input file-> "); scanf("%s",file_name); if((fptr=fopen(file_name,"rb"))==NULL) { printf("file %s does not exist.\n"); exit(1); } nsq=(double)filelength(fileno(fptr)); N=sqrt(nsq); m=(int)(log10((double)N)/log10((double)2.0)); k=1 ; for(i=0;i<m;i++) k<<=1 ; if (k!=N) { printf("Length of file has to be multiples of 2.\n "); exit(1); } Hình 6.12 ảnh dịch chuyển của “MOHSEN.IMG”. printf(" Enter file name for output file.-->"); scanf("%s",file_name); fptro=fopen(file_name,"wb+"); /* 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); /* Allocating memory for twiddle factors. */ n2=(N>>1)-1; wr=(float *)malloc(n2*sizeof(float)); wi=(float *)malloc(n2*sizeof(float)); /*Generating LUT for twiddle factors. */ WTS(wr,wi,N,-1); FFT2D(fptr,fptro,wr,wi,L,N,m,-1); } void FFT2D(FILE *fptr, FILE *fptro, float *wr, float *wi, unsigned int *L, int N, int m, int sign) { /* 2-D FFT Algorithm. */ /* fptr=file pointer to input file. fptro=file pointer to output file. Note: fptr, fptro should be opened in the main routine. They are closed before return to @he main routine. wr[1,wj[I input arrays for twiddle factors, calculated by calling procedure WTS. L[I look-up table for bit reversal. N input array size ( NxN words.) =2 to the power m. sign =-1 for 2-D FFT, =1 for 2-D IFFT. For FFT (sign= 1) the input data is assumed to be real. The result of the FFT has its origin shifted to (N/2,N/2).*/ int N2,i,j,k,kk; long loc,NB; float *xr,*xi,*buff; N2=N<<1; NB=sizeof(float)*N2; /* Allocating memory for FFT arrays ( real and imag.) */ xr=(float *)malloc(N*sizeof(float)); xi=(float *)malloc(N*sizeof(float)); buff=(float *)malloc(NB); /* First stage. */ gotoxy(1,3); printf(" First stage. "); for(j=0;j<N;j++) { gotoxy(1,4); printf("FFT of row %4d",j); if(sign==(int)1) { fread(buff,NB,1,fptr); for(i=0;i<N;i++) { k=L[i] ; kk=i<<1 ; xr[k]=buff[kk]; xi[k]=buff[kk+1]; } } else { for(i=0;i<N;i++) { k=L[i] ; xr[k]=(float)getc(fptr); if(((i+j)%2)!=0) xr[k]=-xr[k]; xi[k]=0.0; } } FFT(xr,xi,wr,wi,m,N); for(i=0;i<N;i++) { k=i<<1; buff[k]=xr[i]; buff[k+1]=xi[1]; } fwrite(buff,NB,1,fptro); } fclose(fptr); /* Transpose. */ gotoxy(1,5); printf(" Transposing of intermediate file. "); rewind(fptro); transpose(fptro,N,m); rewind(fptro); /* Second stage. */ printf("\n Second stage."); for(j=0;j<N;j++) { gotoxy(1,7); printf("FFT of row %4d",j); loc=(long)j*NB; fseek(fptro,loc,SEEK_SET); fread(buff,NB,1,fptro); for(i=0;i<N;i++) { k=L[i] ; kk=i<<1 ; xr[k]=buff[kk]; xi[k]=buff[kk+1]; } FFT(xr,xi,wr,wi,m,N); for(i=0;i<N;i++) { k=i<<1 ; if((sign==(int)1)&&((i+j)%2)==(int)0) { buff[k]=-xr[i]; buff[k+1]=-xi[i]; } else { buff[k]=xr[i]; buff[k+1]=xi[i]; } } if(fseek(fptro,loc,SEEK_SET)!=0) { perror("\n fseek failed.\n"); exit(1) ; } fwrite(buff,NB,1,fptro); } fclose(fptro); } void FFT (float *xr, float *xi, float *wr, float *wi, int m, int N) { /* 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. 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,i,j; float Tr,Ti; ip=1; kk=(N>>1); incr=2 ; for(iter=0; iter<m; iter++) { for(j=0; j<N; 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; } if(iter!=0) { for(k=1; k<ip; k++) { i=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<<=l ; incr<<=l ; } } /********************************************************** This procedure is for obtaining the transpose of a large file of complex numbers stored on secondary memory. We assume that the file is large to the point that it would not fit in active memory. ***********************************************************/ struct COMPLEX { float REAL; float IMAG; } ; void transpose(FILE *fptr, int N, int n) /* Algorithm */ { int N1 , inc ; int iter,i,k; int k1,inc1; int k2,j,k3,k4,NS; struct COMPLEX *buff1,*buff2,tmp; long loc,NT; NS=sizeof(struct COMPLEX); NT=N*NS ; buff1=(struct COMPLEX *)malloc(NT); buff2=(struct COMPLEX *)malloc(NT); N1=N/2 ; inc=1 ; inc1=2 ; for(iter=0;iter<n;iter++) { gotoxy(1,6); printf(" Iteration number %4c.",(iter+1)); k1=0 ; for(k=0;k<N1;k++) { for(i=k1;i<(k1+inc);i++) { loc=NT*(long)(i); if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1); } else fread(buff1,NT,1,fptr); loc=NT*(long)(j+inc); if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1) ; } else fread(buff2,NT,1,fptr); k3=0 ; for(k2=0;k2<N1;k2++) { for(j=k3;j<(k3+inc);j++) { tmp=*(buff1+j+inc); *(buff1+j+inc)=*(buff2+j); *(buff2+j)=tmp; } k3+=inc1 ; } loc=NT*(long)i; if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1); } else fwrite(buff1,NT,1,fptr); loc=NT*(long)(i+inc); if(fseek(fptr,loc,SEEK_SET)!=0) { perror("fseek failed"); exit(1) ; } else fwrite(buff2,NT,1,fptr); } k1+=inc1; } inc*=2 ; inc1*=2 ; N1/=2 ; } } Bài tập 6.5 Cho các mảng 2-D và Phát triển một chương trình C thực hiện Phát triển chương trình C tính tích chập tuần hoàn giữa hai dãy trong miền không gian. Phát triển chương trình C mà sẽ thêm các điểm 0 để mỗi chiều của mảng có độ dài ít nhất là 3 + 3 – 1 = 5 và định dạng tích chập tuần hoàn qua DFT. Dùng chương trình 6.6 để đưa ra tích chập tuần hoàn qua 2-D FFT. 6.6 Hiển thị FFT Nếu FFT của một ảnh trong trường hợp tổng quát là một mảng của các số phức đầy đủ, người ta thường biểu diễn biên độ và pha của tần số của ảnh. Hai yếu tố này biểu diễn tính chất của ảnh. Thông thường biên độ tần số được biểu diễn riêng lẻ và gọi là phổ biên độ. Mặc dù vậy, như chúng ta đã nghiên cứu, pha đóng vai trò quan trọng trong xử lý ảnh, và hợp không hợp lý khi chỉ biểu diễn phổ biên độ của ảnh. Để biểu diễn phổ dưới dạng ảnh, tất cả các việc chúng ta cần phải làm là chia biên độ của FFT thành các giá trị từ 0 đến 255 (cho ảnh 8 bit). Dù thế nào đi chăng nữa thì phổ của ảnh cũng bị suy giảm rất nhanh khi tần số tăng lên. Vì vậy mà vùng tần số cao sẽ trở nên lu mờ khi biểu diễn phổ dưới dạng ảnh. Để giải quyết vấn đề này chúng ta cần xử lý biên độ phổ một chút bằng hàm log. Hàm logarit sẽ sửa độ khuếch đại, và thay thế cho hiển thị phổ |H(u,v)| chúng ta hiển thị: D(u,v) = log10(1+|H(u,v)|) (6.67) Biểu thức này cho ta giá trị zero khi D(u,v) = 0 hay |H(u,v)| = 0 và như vậy D(u,v) luôn luôn có giá trị dương. Một chương trình dùng để chuyển đổi phổ thành dạng ảnh được cho ở chương trình 6.7. Hình 6.13 giới thiệu phổ của ảnh "IKRAM.IMG" trong hình 3.2a sau khi được chuyển đổi dùng biểu thức (6.67). Điểm tần số (0,0) nằm ở trung tâm màn hình. Chú ý phổ ảnh giảm xuống rất nhanh chóng khi tần số tăng lên. Chương trình 6.7 " DISP_FFT.C" Chương trình dùng để đưa ra một file chứa phổ tần số trong dạng ảnh có thể hiển thị được. /************************ * Program developed by: * * M.A.Sid-Ahmed. * * ver.1.0 1992. * *************************/ /**************************************************** Program for calculating the magnitude of the 2-D FFT given a file containing the complex values of the FFT of an image. The result is placed in a form suitable for display in image form and stored in an external file. The mapping function D(u,v)=log10(1+ |(F(u,v)|) is used. ******************************************************/ #include #include #include #include #include #include void main() { int i,j,k,N,NB1,NB2; FILE *fptri, *fptro,*fptrt; double nsq; float max,min,xr,xi,scale; float *buffi,*bufft; char *buffo; char file_name[14]; Hình 6.13 Phổ của "IKRAM.IMG" clrscr(); printf(""Enter name of file containing FFT data ---> "); scanf("%s",file_name); fptri=fopen(file_name,"rb"); printf("Enter name of file for storing magnitude data -> "); scanf("%s",file_name); fptro=fopen(file_name,wb"); fptrt=fopen("temp.img","wb+"); nsq=(double)(filelength(fileno(fptri))/(2*sizeof(float))); N=(int)(sprt(nsq)); max=0.0; min=1.0e9; NB1=(N<<1)*sizeof(float); NB2=NB1>>1; buffi=(float *)malloc(NB1); bufft=(float *)malloc(NB2); buffo(char *)malloc(N*sizeof(char)); for(i=0;i<N;i++) { fread(buffi,NB1,1,fptri); for(j=0;j<N;j++) { k=j<<1; xr=buffi[k]; xi=buffi[k+1]; bufft[j]=(float)sqrt((double)(xr*xr+xi*xi)); bufft[j]=(float)log10((double)(1+bufft[j])); if(bufft[j]>max) max=bufft[j]; if(bufft[j]<min min=bufft[j]; } fwrite(bufft,NB2,1,fptrt); } fclose(fptri); fseek(fptrt,0,SEEK_SET); scale=(float)255.0/(max-min); for(i=0;i<N;i++) { fread(bufft,NB2,1,fptrt); for(j=0;j<N;j++) buffo[j]=(char)((bufft[j]-min)*scale); fwrite(buffo,N,1,fptro); } fclose(fptro); fclose(fptrt); remove("temp.img"); } Bài tập 6.6 Tính 2D_FFT của " IKRAM.IMG" và hiển thị |H(u,v)| thay thế cho hiện thị D(u,v). Chú ý dến sự suy giảm phổ ảnh và so sánh với trường hợp hiển thị ảnh rút ra bởi chương trình 6.7. Bài tập 6.7 Lập một chương trình 2-D FFT theo các bước sau: 1. Thuật toán phân chia tần số. 2. Thuật toán giảm lược đầu vào. 3. Thật toán giảm lược đầu ra. 4. Dùng thuật toán giảm lược đầu ra thiết kế một bộ lọc 2-D FIR thông thấp vói D0 = 0.3p, kích thước 11 ´ 11. So sánh ví dụ 2.5 trong chương 2. 6.7 Bộ lọc hai chiều dùng FFT Nếu dùng tích chập để chuyển hàng loạt các phần tử từ miền không gian sang miền tần số ta nên áp dụng FFT. Phép biến đổi này yêu cầu 2. (N2/2). log2N phép nhân phức và 2. N2. log2N phép cộng phức để thu được 2-D FFT, N2 phép nhân phức trong miền tần số giữa FFT của điểm ảnh và các đáp ứng tần số cuả bộ lọc, 2 . (N2/2) . log2N phép nhân phức cho IFFT. Mặt khác, một bộ lọc 2-D FIR có kích thước (2m + 1) ´ (2m + 1) đòi hỏi (2m + 1)2 N2 phép nhân để thu được ảnh trực tiếp trong miền không gian. Xem xét một ảnh có kích thước 512 ´ 512 điểm. FFT yêu cầu: ằ 20 triệu phép nhân. Để đưa ra tính toán này chúng ta coi rằng một phép nhân phức thì bằng 4 phép nhân thông thường, và bộ lọc có pha zero. Phương pháp không gian áp dụng cho một bộ lọc có kích thước 7 ´ 7 yêu cầu 7 ´ 7 ´ 5122 ằ 13 triệu phép nhân. Nếu kích thước bộ lọc tăng lên thì phương pháp phân chia miền tần số có thể áp dụng. Một bộ lọc có kích thước 11 ´ 11 yêu cầu khoảng 30 triệu phép nhân sẽ chỉ cần khoảng 19 triệu phép nhân khi áp dụng phương pháp phân chia miền tần số. Hai phương pháp này sẽ có cùng một số phép nhân nếu Cho một ảnh có kích thước 512 ´ 512 (2m + 1) ằ 9, dễ chứng minh là nếu kích thước bộ lọc nhỏ hơn 9 thì ta có thể phương pháp phân chia không gian. Tuy nhiên, cần chú ý phương pháp phân chia tần số cũng yêu cầu ít thời gian xử lý hơn do số lần truy nhập đĩa giảm xuống. Ưu điểm này được tăng lên khi kích thước của bộ lọc lớn hơn 9 ´ 9. Ưu điểm này sẽ không còn nữa khi xét đến lỗi wraparound. Để tránh lỗi này ta phải tăng gấp bốn lần kích thước của ảnh. Cho một ảnh có kích thước 512 ´ 512 ta cần phải tăng lên 1024 ´ 1024. Để tránh các phép tính toán quá lớn khi chú ý rằng h(n1, n2) của một bộ lọc khi rút ra IFFT sẽ tăng lên rất nhanh khi n1, n2 tăng lên. Tính chất này càng nổi bật khi mở rộng Fourier chỉ chèn các giá trị zero vào các giá trị cuối của bộ lọc từ . Cần nhắc lại là cả đáp ứng tấn số và đáp ứng xung được xem xét khi làm việc với DFT. Thuộc tính là h(n1, n2) tăng lên một cách nhanh chóng được xem xét khi lựa chọn phương án lọc. Không phụ thuộc vào kích thước của ảnh, đưa ra phép nhân giứa đáp ứng tần số của ảnh và đáp ứng tần số của bộ lọc, và chúng ta chú ý rằng lỗi wrapapound chỉ xuất hiện ở miền nhỏ nằm ở đường bao của ảnh và trong phần lớn trường hợp lỗi này có thể bỏ qua. Phương pháp tần số có thể thực hiện qua các bước sau: Rút ra 2-D FFT của một ảnh Nhân I(k1, k2) với đặc tuyến của bộ lọc, chú ý là đáp ứng tần số có gốc toạ độ nằm tại (N/2, N/2). Cho ví dụ một bộ lọc thông cao Butterworth có đặc tuyến như sau: Dùng Đáp ứng tần số của ảnh lọc có thể rút ra từ ảnh đã lọc có thể rút ra từ : ở đây  có nghĩa là phần thực của phần nằm trong hai dấu ngoặc. Bài tập 6.8 Viết một chương trình lọc 2-D trong mặt phẳng tần số. Kiểm tra chương trình dùng cùng một đặc tuyến tần số dùng thiết kế bộ lọc FIR lọc ảnh trong hình 2.3 (chương 3) và so sánh kết quả. 6.8 Vector biến đổi Fourier Qua chiến lược chia để trị ta đạt được hiệu suất tính toán trên máy tính của giải thuật 1-D FFT. Thuật toán FFT vector 2-D sau đâylà cùng một chiến lược. Giải thuật DFT 2-D được xen kẽ với những giải thuật DFT 2-D nhỏ hơn, cuối cùng chỉ DFT 2-D của phần tử đơn được tính. Chúng ta sẽ kiểm tra vector FFT. Chúng ta có ở đây Vì vậy có thể viết lại như sau: Tương tự với trường hợp1-D chúng ta có thể viết (6.68a) (6.68b) (6.68c) (6.68d) Dạng công thức cuối cùng được biết đến như một bướm cơ số (2 ´ 2). Mỗi bướm cần đến ba phép nhân và tám phép cộng, để tính toán một bước của giải thuật FFT đòi hỏi N2/4 bướm (việc kiểm chứng được giành cho độc giả coi như là một bài tập). Số bước cần thiết để thực hiện giải thuật FFT 2-D là log2N, vì vậy số phép nhân cần thực hiện là Với phương pháp này số hàng-cột cần thiết sẽ nhỏ hơn 25 phần trăm so với số hàng-cột được mô tả trước đây. Tuy nhiên phương pháp này sẽ chỉ có hiệu quả nếu có đủ bộ nhớ hoạt động lưu giữ N ´ N số phức. Bài tập 6.9 1. Phát triển thuật toán và chương trình C cho giải thuật phân chia thời gian 2-D FFT cơ số vector 2. Biến đổi công thức để chia tần số vector FFT 2-D. 3. Phát triển thuật toán giảm lược đầu vào và đầu ra, viết chương trình C cho 2-D FFT cơ số vector.