Tìm hiểu về biến đổi Fourier FFT/IFFT và phương pháp tính toán

byVi-et Spaces 0 تعليقات
Trong các bài tóan phân tích tín hiệu thực tế, tín hiệu cần phân tích có khỏang thời gian xác định rất dài hoặc vô hạn, khi đó phân tích dùng chuỗi Fourier (Fourier series) trở nên không thích hợp. Trong trường hợp này, biến đổi Fourier (Fourier transform (FT)) và biến đổi ngược của nó (inverse Fourier transform (IFT)) được sử dụng. Biến đổi Fourier được áp dụng thành công trong các lĩnh vực khoa học kỹ thuật ở đó khái niệm tần số và miền tần số được sử dụng. Biến đổi Fourier được mở rộng từ chuỗi Fourier bằng cách đặt khỏang thời gian xác định của chuỗi Fourier vô hạn, hoặc biến đổi Fourier có thể định nghĩa độc lập và chuỗi Fourier là trường hợp đặc biệt của nó. 


Definition of the Fourier Transform

Biến đổi Fourier(FT) của hàm thực (hoặc phức) của biến thực t được định nghĩa bởi

X(ω)=∫∞−∞x(t)e−jωtⅆt   (1)

cho hàm phức của biến thực ω biểu diễn tần số. Biến đổi Fourier ngược (IFT) được cho bởi (IFT)
x(t)=12π∫∞−∞X(ω)ejωtⅆω.   (2)

Vì các cận vô hạn của hai tích phân nên vấn đề hội tụ trở nên quan trọng, do đó giới hạn khả năng biểu diễn một số lọai tín hiệu. Sử dụng hàm delta (phân bố) trong miền thời gian và tần số cho phép biểu diễn nhiều lọai tín hiệu. [1]. 
Examples of the Fourier Transform

Xác định các phép biến đổi cơ bản và dùng các tính chất của phép biến đổi Fourier cho phép khảo sát nhiều lọai tín hiệu (ví dụ: điều chế, lấy mẫu) một cách dễ dàng.
  • Nếu x(t)=δ(t) thì X(ω)=1

  • Nếu x(t)=1 thì X(ω)=2πδ(ω)

  • Nếu x(t) là chuỗi vô hạn của các hàm delta đặt cách nhau khỏang T , x(t)=∑∞n=−∞δ(t−nT), biếb đổi của nó cũng là chuỗi vô hạn các hàm delta có trọng số 2π/T đặt cách nhau khỏang 2π/T ,X(ω)=2π∑∞k=−∞δ(ω−2πk/T).

  • Các ví dụ khác có thể tham khảo trong [1][2].
Lưu ý rằng, biến đổi Fourier của hàm liên tục thời gian tạo nên hàm liên tục tần số. Nếu hàm delta được dùng, biến đổi Fourier của hàm tuần hòan sẽ là chuỗi vô hạn của hàm delta với trọng số là hệ số của chuỗi Fourier.

Trong xử lý ảnh số, có những công cụ nhất định giữ vai trò trung tâm. Chúng bao gồm các công cụ toán như tích chập, phân tích Fourier, và những công cụ mô tả và xử lý thống kê như mã dây chuyền và mã chạy. Chúng tôi sẽ trình bày những công cụ này mà không đề cập đến bất kì đẫn giải nào. Phần dẫn giải sẽ xuất hiện ở những mục tiếp sau.

3.1 Tích chập

Có một vài cách kí hiệu để chỉ tích chập của hai tín hiệu (nhiều chiều) để tạo ra một tín hiệu đầu ra. Những kí hiệu thường dùng nhất bao gồm:

c = a ⊗ b = a∗b                  (1)

Ta sẽ dùng dạng thứ nhất, c = a ⊗ b, trong những định nghĩa chính thức sau.

Trong không gian 2 chiều liên tục:

c(x, y) = a(x, y) ⊗ b(x, y) = ∫–∞∞∫–∞∞a(χ, ζ) b(x – χ, y – ζ) dχdζ           (2)

Trong không gian 2 chiều rời rạc:

c[m, n] = a[m, n] ⊗ b[m, n] = ∑j=–∞+∞∑k=–∞+∞a[j, k] b[m – j, n – k]              (3)

3.2 Những đặc tính của tích chập

Có một số đặc tính quan trọng gắn với tích chập.
  • Tích chập có tính giao hoán.
c = a ⊗ b = b ⊗ a                  (4)
  • Tích chập có tính kết hợp.
c = a ⊗ (b ⊗ d) = (a ⊗ b) ⊗ d = a ⊗ b ⊗ d                (5)
  • Tích chập có tính phân phối.
c = a ⊗ (b + d) = (a ⊗ b) + (a ⊗ d)                      (6)

trong đó a, b, c, và d là các ảnh, kiểu liên tục hoặc rời rạc.

3.3 Biến đổi Fourier

Phép biến đổi Fourier cho ta một cách biểu diễn khác cho tín hiệu, cụ thể là biểu diễn dưới dạng một tổng trọng số của các số mũ phức. Từ công thức Euler:

ejq = cosq + j sinq                        (7)

trong đó j2 = –1, ta có thể nói rằng phép biến đổi Fourier cho ta một cách biểu diễn của tín hiệu (2 chiều) dưới dạng một tổng trọng số của các hàm sin và cosin. Những công thức định nghĩa cho các phép biến đổi Fourier thuận và Fourier nghịch như sau. Cho trước một ảnh a và dạng biến đổi Fourier của nó,A, thì phép biến đổi thuận đi từ miền không gian (liên tục hoặc rời rạc) tới miền tần số, vốn luôn có tính liên tục.

Thuận:                A = \mathfrak{<em>F</em>} {a}                               (8)

Phép biến đổi Fourier nghịch đi từ miền tần số ngược về miền không gian.

Nghịch:               a = F{A}                              (9)

Phép biến đổi Fourier là phép toán duy nhất và không thể đảo ngược, sao cho:

a = F-1{F{a}}   và   A = F{F-1{A}}                (10)

Công thức cụ thể để biến đổi thuận và nghịch giữa miền không gian và miền tần số được cho sau đây:

Trong không gian 2 chiều liên tục:

Thuận –                   A(u, v) = ∫–∞+∞∫–∞+∞a(x, y) e–j(ux+vy)dxdy                           (11)

Nghịch –                 a(x, y) = (1/4π2) ∫–∞+∞∫–∞+∞A(u, v) e+j(ux+vy)dudv        (12)

Trong không gian 2 chiều rời rạc: 

Thuận –                  A(Ω, Ψ) = ∑m=–∞+∞∑n=–∞+∞a[m, n]e–j(Ωm+Ψn) 

Nghịch –               a[m, n] = (1/4π2) ∫–ππ∫–ππA(Ω, Ψ)e+j(Ωm+Ψn)dΩdΨ

3.4 Những đặc tính của biến đổi Fourier

Có một loạt những đặc tính gắn với biến đổi Fourier thuận và nghịch. Dưới đây là một số đặc tính quan trọng nhất đối với xử lý ảnh:
  • Nói chung, phép biến đổi Fourier là một hàm của các biến tần số thực. Theo đó phép biến đổi có thể được viết dưới dạng độ lớn và pha của nó.
A(u, v) = |A(u, v)|e jϕ(u, v)              A(Ω, Ψ) = |A(Ω, Ψ)|e jϕ(Ω, Ψ)                 (15)
  • Một tín hiệu 2 chiều cũng có thể có dạng phức và được viết dưới dạng độ lớn và pha của nó
a(x, y) = |a(x, y)|e jϑ(x, y)               a[m, n]=|a[m, n]|e jϑ[m, n]                      (16)

  • Nếu tín hiệu hai chiều có dạng thực, thì biến đổi Fourier sẽ có những cặp giá trị đối xứng nhất định:
A(u, v) = A∗(–u, –v)                        A(Ω,Ψ) = A∗(–Ω, –Ψ)                            (17)

Kí hiệu (∗) để chỉ liên hợp phức. Với các tín hiệu thực, PT (17) sẽ trực tiếp dẫn tới: 

|A(u, v)| = |A(–u, –v)|          ϕ(u, v)= –ϕ(–u, –v)
|A(Ω, Ψ)| = |A(–Ω, –Ψ)|     ϕ(Ω, Ψ)= –ϕ(–Ω, –Ψ)                     (18)

  • Nếu tín hiệu hai chiều có tính thực và chẵn, thì biến đổi Fourier của nó cũng là thực và chẵn.
A(u, v) = A(–u, –v)                A(Ω, Ψ) = A(–Ω, –Ψ)                      (19)
  • Các phép biến đổi Fourier thuận và nghịch là các phép toán tuyến tính.
F{w1a + w2b} = F{w1a} + F{w2b} = w1A + w2B
F–1{w1a + w2b} = F–1{w1a} + F–1{w2b} = w1A + w2B               (20)

trong đó a và b là các tín hiệu (ảnh) 2 chiều còn w1 và w2 là các hằng số phức tùy ý.
  • Phép biến đổi Fourier trong không gian rời rạc, A(Ω,Ψ), có tính chu kỳ trên cả Ω và Ψ. Cả hai chu kì đều là 2π.
A(Ω + 2πj, Ψ + 2πk) = A(Ω, Ψ)              j,k nguyên                              (21)
  • Năng lượng, E, trong một tín hiệu có thể được đo trên miền không gian hoặc trên miền thời gian. Với một ín hiệu có nặng lượng hữu hạn:Định lý Parseval (không gian 2 chiều liên tục):
E = ∫–∞+∞∫–∞+∞|a(x,y)|2dxdy = (1/4π2) ∫–∞+∞∫–∞+∞|A(u,v)|2dudv                 (22)

Định lý Parseval (không gian 2 chiều rời rạc):

E = ∑m=–∞+∞ ∑n=–∞+∞|a[m, n]|2 = (1/4π2) ∫–π+π∫–π+π|A(Ω, Ψ)|2dΩ dΨ      (23)

Không nên nhầm “năng lượng tín hiệu” này với nặng lượng vật lý ở trong hiện tượng phát ra tín hiệu. Chẳng hạn, nếu giá trị a[m,n] biểu thị số các photon, thì năng lượng vật lý sẽ tỉ lệ với biên độ a, chứ không phải là với bình phương của biên độ. Nói chung, trường hợp này được xét với việc ghi hình video.
  • Cho trước ba tín hiệu 3 chiều a, b, và c cùng các dạng biến đổi Fourier tương ứng của chúng, A, B, C:
c = a ⊗ b ←F→ C = A • B

c = a • b ←F→ C = (1/4π2) A ⊗ B                   (24)

Diễn đạt bằng lời thì tích chập trong miền không gian tương đương với phép nhân trong miền tần số (miền Fourier) và ngược lại. Đây là một kết quả trung tâm; nó không chỉ cung cấp một phương pháp thực hiện tích chập mà còn làm rõ cách mà hai tín hiệu tương tác với nhau—bằng cách chập—để hình thành một tín hiệu thứ ba. Ta sẽ còn sử dụng nhiều đến kết quả này.

  • Nếu một tín hiệu 2 chiều a(x,y) được phóng đại theo tọa độ không gian của nó thì:
    Nếu    a(x, y) → a(Mx•x, My•y)
    Thì    A(u, v) → A(u/Mx, v/My)/|Mx•My|                                   (25)

  • Nếu một tín hiệu 2 chiều a(x,y) có phổ Fourier A(u,v) thì:
    A(u=0, v=0) = ∫–∞+∞∫–∞+∞a(x,y)dxdy
    a(x=0,y=0) = (1/4π2) ∫–∞+∞∫–∞+∞A(u,v)dxdy                      (26)

  • Nếu một tín hiệu 2 chiều a(x,y) có phổ Fourier A(u,v) thì:
    ∂a(x,y)/∂x ←F→ juA(u,v)             ∂a(x,y)/∂y ←F→ jvA(u,v)
    ∂2a(x,y)/∂x2 ←F→ –u2A(u,v)     ∂2a(x,y)/∂y2 ←F→ –v2A(u,v)            (27)

TẦM QUAN TRỌNG CỦA PHA VÀ BIÊN ĐỘ

Phương trình (15) cho thấy rằng biến đổi Fourier của một hình ảnh có thể sẽ phức tạp. Điều đó được minh họa ở Hình 4a-c sau đây. Hình 4a cho thấy ảnh gốc a[m,n]. Hình 4b là biên đồ được giãn theo tỉ lệ log(|A(Ω,Ψ)|), còn Hình 4c là pha ϕ(Ω,Ψ).


Hình 4 

Cả hai hàm biên độ và pha đều cần có để tái tạo hoàn toàn được ảnh từ dạng biến đổi Fourier của nó. Hình 5a cho thấy điều gì sẽ sảy ra nếu Hình 4a được tái tạo chỉ dựa trên thông tin về biên độ, còn Hình 5b cho thấy điều gì sẽ xảy a nếu Hình 4a được tái tạo chỉ dựa trên cơ sở thông tin về pha.

Hình 5 

Như vậy, chỉ với thông tin về biên độ hoặc pha không thôi thì vẫn chưa đủ để tái tạo lại hình. Hình chỉ tạo với biên độ (Hình 5a) không thể nhận ra được và có vấn đề về khoảng sáng tối (dynamic range). Hình chỉ tạo với pha (Hình 5b) thì chỉ đủ để nhận ra, nghĩa là có chất lượng xuống cấp nghiêm trọng.

TÍN HIỆU ĐỐI XỨNG VÒNG

Một tín hiệu 2 chiều a(x,y) bất kì luôn có thể được biểu diễn trong hệ tọa độ dưới dạng a(r,θ). Khi tín hiệu 2 chiều này có tính đối xứng vòng, điều này nghĩa là:

a(x,y) = a(r, θ) = a(r)                         (28)

trong đó r2 = x2 + y2 và tanθ = y/x. Vì một loạt các hệ vật lý như thấu kính thể hiện tính đối xứng vòng, nên sẽ rất có ích nếu ta tính được dạng biến đổi Fourier phù hợp.

Dạng biến đổi Fourier A(u,v) có thể được viết trong hệ tọa độ cực là A(q,ξ) và theo đó, đối với một tín hiệu có tính đối xứng vòng, được viết dưới dạng biến đổi Hankel:

A(u,v) = F{a(x,y)} = 2π∫0∞a(r)J0(rq)rdr = A(q)         (29)

trong đó q2 = u2 + v2, tanξ = v/u còn J0(•) là hàm Bessel dạng thứ nhất với bậc bằng 0.

Dạng biến đổi Hankel ngược được cho bởi:

a(r) = 1/2π ∫0∞A(q)J0(rq)qdq                                            (30)

Dạng biến đổi Fourier của một tín hiệu 2 chiều có tính đối xứng vòng là một hàm của biến duy nhất là tần số dài, q. Sự phụ thuộc của tần số góc ξ đã không còn. Hơn nữa, nếu a(x,y) = a(r) là hàm số thực thì nó sẽ mặc nhiên là hàm chẵn do tính đối xứng vòng. Theo Phương trình (19), A(q) cũng sẽ là hàm số thực và chẵn.



3.4.3 VÍ DỤ VỀ CÁC TÍN HIỆU 2 CHIỀU CÙNG DẠNG BIẾN ĐỔI CỦA CHÚNG


Bảng 4 cho ta thấy một số tín hiệu cơ bản và hữu ích cùng với các biến đổi Fourier 2 chiều của chúng. Khi dụng các dạng cho trong bảng trong phần còn lại của chương này, ta sẽ gọi biểu thức trong miền không gian là hàm rải điểm (point spread function, PSF) hay phản hồi xung 2 chiều (2D impulse response) và dạng biến đổi Fourier của nó là hàm biến đổi quang học (optical transfer function, OTF), hay đơn giản là hàm biến đổi. Hai tín hiệu tiêu chuẩn được dùng trong bảng này là u(•), và hàm bậc đơn vị, và J1(•), hàm Bessel dạng thứ nhất. Các tín hiệu đối xứng vòng được coi như là các hàm theo rnhư trong PT (28).

image image image image image image image image

Bảng 4. Các ảnh 2 chiều cùng dạng biến đổi Fourier của chúng

Phép biến đổi Fourier nhanh


Phương pháp được phác họa trong Mục [ĐK biên Dirichlet ở Phần 1] để giải phương trình Poisson trong 2 chiều với điều kiện biên Dirichlet đơn giản theo phương y đòi hỏi ta phải thực hiện rất nhiều lượt biến đổi Fourier cho hàm sin:

FjS = 2 / J ∑ k = 1J − 1fk sin(j k π / J)  (ffta)

với j = 0, J, và các biến đổi Fourier ngược cho hàm sin:

fj = ∑ k = 1J − 1FkS sin(j k π / J).   (fftb)

Ở đây, fj là giá trị của f(y) tại yj = j L / J. Như vậy, PT ([ffta]) giống như ([ffta1]) và ([ffta2]), còn PT ([fftb]) có thể được dùng để tái lập u(xi, yj) từ cácUi, j. Tương tự, phương pháp được vạch ra trong Mục [ĐK biên Neumann ở Phần 1] để giải hệ phương trình Poisson trong 2 chiều với các điều kiện biên Neumann đơn giản theo phương y đòi hỏi chúng ta phải thực hiện rất nhiều lượt biến đổi Fourier cho hàm cosin:

\displaystyle{ F_j^C = \frac{f_0}{J} + \frac{2}{J}\sum_{k=1}^{J-1} f_k\,\cos(j\,k\,\pi/J) + \frac{(-1)^j\,f_J}{J} }

với j = 0, J, và các biến đổi Fourier ngược cho hàm cosin:

fj = ∑ k = 0JFkC cos(j k π / J). 

Thật không may là để thực hiện trực tiếp các biến đổi như vậy, cần đến O(J2) phép tính số học, nghĩa là chúng cực kỳ tiêu tốn tài nguyên CPU. Tuy nhiên, có một thuật toán khéo léo dành riêng cho việc biến đổi Fourier mà chỉ tốn O(J lnJ) phép tính số học [vốn nhỏ hơn O(J2) đáng kể khi J lớn]. Thuật toán này có tên là biến đổi Fourier nhanh hay FFT.1

Chi tiết của thuật toán FFT vượt ra ngoài phạm vi khóa học này. Nói nôm na, thuật toán thực hiện việc chuyển đổi theo các giai đoạn nối tiếp dùng 2, 4, 8, 16, v.v. điểm nút. Ở khóa học này, ta sẽ dùng thư viện FFT thông dụng nhất có thể kiếm được trên mạng, thư viện fftw,2 để thực hiện tất cả các phép biến đổi Fourier cho hàm sin và cosin. Không may là thư viện fftw không trực tiếp thực hiện tính toán biến đổi Fourier các hàm sin và cosin.3 Thay vì vậy, nó tính biến đổi Fourier phức:

Fj = 1 / 2J ∑ k = 02J − 1fk exp( − i j k π / J)

với j = 0, 2J − 1, và các biến đổi Fourier ngược phức:

fj = ∑ k = 02J − 1Fk exp( i j k π / J). 

Lưu ý rằng fj và Fj có tính tuần hoàn theo j với chu kì 2 J. Cũng lưu ý rằng bộ số liệu gắn liền với phép biến đổi Fourier phức bao gồm số phần tử nhiều gấp hai lần số phần tử trong biến đổi sin hoặc cosin. Tuy nhiên, ta có thể dễ dàng chuyển từ phép biến đổi sin hoặc cosin sang biến đổi phức bằng cách kéo dàibộ số liệu. Vì vậy, đối với biến đổi sin ta viết:

\begin{array}{ccc} f_{2J-j} &=& - f_j,\\[0.5ex] F_{2J-j} &=& -F_j, \end{array}

với j = 1, J − 1, trong trường hợp này

FjS = 2 i Fj

Tương tự, với biến đổi cosin ta viết:

\begin{array}{ccc} f_{2J-j} &=& f_j,\\[0.5ex] F_{2J-j} &=&F_j, \end{array}

với j = 1, J − 1, trong trường hợp này

FjC = 2 Fj. 

Dưới đây là một loạt các đoạn chương trình để bọc các hàm trong thư viện fftw để thực hiện biến đổi Fourier cho các hàm sin và cosin.
// FFT.cpp 

// Set of functions to calculate Fourier-cosine and -sine transforms 
// of real data using fftw Fast-Fourier-Transform library. 
// Input/ouput arrays are assumed to be of extent J+1. 
// Uses version 2 of fftw library (incompatible with vs 3). 

#include <fftw.h> 
#include <blitz/array.h> 

using namespace blitz; 

// Calculates Fourier-cosine transform of array f in array F 
void fft_forward_cos (Array<double,1> f, Array<double,1>& F) 
{ 
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (ff[0]) = f(0); c_im (ff[0]) = 0.; 
  c_re (ff[J]) = f(J); c_im (ff[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.; 
      c_re (ff[2*J-j]) = f(j); c_im (ff[2*J-j]) = 0.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_one (p, ff, FF); 
  fftw_destroy_plan (p);  

  // Unload data 
  F(0) = c_re (FF[0]); F(J) = c_re (FF[J]);  
  for (int j = 1; j < J; j++) 
    { 
      F(j) = 2. * c_re (FF[j]); 
    } 

  // Normalize data 
  F /= 2. * double (J); 
} 

// Calculates inverse Fourier-cosine transform of array F in array f 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f) 
{     
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (FF[0]) = F(0); c_im (FF[0]) = 0.; 
  c_re (FF[J]) = F(J); c_im (FF[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (FF[j]) = F(j) / 2.; c_im (FF[j]) = 0.; 
      FF[2*J-j] = FF[j]; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE); 
  fftw_one (p, FF, ff); 
  fftw_destroy_plan (p);  

  // Unload data  
  f(0) = c_re (ff[0]); f(J) = c_re (ff[J]);  
  for (int j = 1; j < J; j++) 
    { 
      f(j) = c_re (ff[j]); 
    } 
} 

// Calculates Fourier-sine transform of array f in array F 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F) 
{   
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data  
  c_re (ff[0]) = 0.; c_im (ff[0]) = 0.;  
  c_re (ff[J]) = 0.; c_im (ff[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.; 
      c_re (ff[2*J-j]) = - f(j); c_im (ff[2*J-j]) = 0.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_one (p, ff, FF); 
  fftw_destroy_plan (p);  

  // Unload data 
  F(0) = 0.; F(J) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      F(j) = - 2. * c_im (FF[j]); 
    }   

  // Normalize data 
  F /= 2. * double (J); 
} 

// Calculates inverse Fourier-sine transform of array F in array f 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f) 
{  
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (FF[0]) = 0.; c_im (FF[0]) = 0.;  
  c_re (FF[J]) = 0.; c_im (FF[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (FF[j]) = 0.; c_im (FF[j]) = - F(j) / 2.; 
      c_re (FF[2*J-j]) = 0.; c_im (FF[2*J-j]) = F(j) / 2.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE); 
  fftw_one (p, FF, ff); 
  fftw_destroy_plan (p);  

  // Unload data 
  f(0) = 0.; f(J) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      f(j) = c_re (ff[j]); 
    } 
}

Ví dụ tính toán 2 chiều cho bài toán Poisson 2 chiều


Dưới đây là một ví dụ đoạn chương trình giải bài toán Poisson 2 chiều có dùng các đoạn chương trình làm nhiệm vụ nghịch đảo ma trận ba đường chéo và bao bọc các hàm FFT nêu trên, cũng như thư viện Blitz++.
// Poisson2d.cpp 

// Function to solve Poisson's equation in 2-d: 

//  d^2 u / dx^2 + d^2 u / dy^2 = v  for  xl <= x <= xh  and  0 <= y <= L 

//  alphaL u + betaL du/dx = gammaL(y)  at x=xl 

//  alphaH u + betaH du/dx = gammaH(y)  at x=xh 

// In y-direction, either simple Dirichlet boundary conditions: 

//  u(x,0) = u(x,L) = 0 

// or simple Neumann boundary conditions: 

//  du/dy(x,0) = du/dy(x,L) = 0 

// Matrices u and v assumed to be of extent N+2, J+1. 
// Arrays gammaL, gammaH assumed to be of extent J+1. 

// Now, (i,j)th elements of matrices correspond to 

//  x_i = xl + i * dx    i=0,N+1 

//  y_j = j * L / J      j=0,J 

// Here, dx = (xh - xl) / (N+1) is grid spacing in x-direction. 

// Now, kappa = pi * dx / L 

// Finally, Neumann=0/1 selects Dirichlet/Neumann bcs in y-direction. 

#include <blitz/array.h> 

using namespace blitz; 

void fft_forward_cos (Array<double,1> f, Array<double,1>& F); 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f); 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F); 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f); 
void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,  
                  Array<double,1> w, Array<double,1>& u); 

void Poisson2D (Array<double,2>& u, Array<double,2> v,  
                double alphaL, double betaL, Array<double,1> gammaL,  
                double alphaH, double betaH, Array<double,1> gammaH, 
                double dx, double kappa, int Neumann) 
{ 
  // Find N and J. Declare local arrays. 
  int N = u.extent(0) - 2; 
  int J = u.extent(1) - 1; 
  Array<double,2> V(N+2, J+1), U(N+2, J+1); 
  Array<double,1> GammaL(J+1), GammaH(J+1); 

  // Fourier transform boundary conditions 
  if (Neumann) 
    { 
      fft_forward_cos (gammaL, GammaL); 
      fft_forward_cos (gammaH, GammaH); 
    } 
  else 
    { 
      fft_forward_sin (gammaL, GammaL); 
      fft_forward_sin (gammaH, GammaH); 
    } 

  // Fourier transform source term 
  for (int i = 1; i <= N; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = v(i, j); 

      if (Neumann) 
        fft_forward_cos (In, Out); 
      else 
        fft_forward_sin (In, Out); 

      for (int j = 0; j <= J; j++) V(i, j) = Out(j); 
    } 

  // Solve tridiagonal matrix equations 
  if (Neumann) 
    { 
      for (int j = 0; j <= J; j++) 
        {   
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), uu(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = 1.; 
          for (int i = 1; i <= N; i++) 
            b(i) = -2. - double (j * j) * kappa * kappa; 
          b(1) -= betaL / (alphaL * dx - betaL); 
          b(N) += betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = 1.; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
            w(i) = V(i, j) * dx * dx; 
          w(1) -= GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) -= GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, uu); 
          for (int i = 1; i <= N; i++) U(i, j) = uu(i); 
        } 
    } 
  else 
    { 
      for (int j = 1; j < J; j++) 
        {  
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), uu(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = 1.; 
          for (int i = 1; i <= N; i++) 
            b(i) = -2. - double (j * j) * kappa * kappa; 
          b(1) -= betaL / (alphaL * dx - betaL); 
          b(N) += betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = 1.; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
            w(i) = V(i, j) * dx * dx; 
          w(1) -= GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) -= GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, uu); 
          for (int i = 1; i <= N; i++) U(i, j) = uu(i); 
        } 

      for (int i = 1; i <= N  ; i++)  
        { 
          U(i, 0) = 0.; U(i, J) = 0.; 
        } 
    } 

  // Reconstruct solution via inverse Fourier transform 
  for (int i = 1; i <= N; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = U(i, j); 

      if (Neumann) 
        fft_backward_cos (In, Out); 
      else 
        fft_backward_sin (In, Out); 

      for (int j = 0; j <= J; j++) u(i, j) = Out(j); 
    }   // Calculate i=0 and i=N+1 values   
  for (int j = 0; j <= J; j++) 
    { 
      u(0, j) = (gammaL(j) * dx - betaL * u(1, j)) / 
        (alphaL * dx - betaL); 
      u(N+1, j) = (gammaH(j) * dx + betaH * u(N, j)) / 
        (alphaH * dx + betaH); 
    } 
}

Một ví dụ giải phương trình Poisson 2 chiều


Bây giờ ta hãy dùng kĩ thuật nêu trên để giải phương trình Poisson trong không gian hai chiều. Giả sử rằng số hạng nguồn là

v(x, y) = 6 x y (1 − y) − 2 x3

với 0 ≤ x ≤ 1 và 0 ≤ y ≤ 1. Các điều kiện biên tại x = 0 là αl = 1, βl = 0, và γl = 0 [xem PT ([e527a])], trong đó các điều kiện biên tại x = 1 là αh = 1,βh = 0, và γh = y (1 − y) [xem PT ([e528a])]. Các điều kiện biên Dirichlet đơn giản u(x, 0) = u(x, 1) = 0 được áp dụng tại y = 0 và y = 1. Dĩ nhiên là bài toán này có thể giải được theo cách giải tích để cho ta

u(x, y) = y (1 − y) x3. 

Các Hình [p2da] và [p2db] so sánh giữa các nghiệm giải tích và sai phân hữu hạn cho N = J = 64. Có thể thấy rằng nghiệm sai phân hữu hạn khá trùng khớp với nghiệm giải tích.

p2daHình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trụcx tại y = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2da]

p2db

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trụcy tại x = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2db]

Xét ví dụ thứ hai, giả sử rằng số hạng nguồn là

v(x, y) =  − 2 (2 y3 − 3 y2 + 1) + 6 (1 − x2) (2 y − 1)

với 0 ≤ x ≤ 1 và 0 ≤ y ≤ 1. Các điều kiện biên tại x = 0 là αl = 1, βl = 0, và γl = 2 y3 − 3 y2 + 1 [xem PT ([e527a])], còn các điều kiện biên tại x = 1 làαh = 1, βh = 0, và γh = 0 [xem PT ([e528a])]. Các điều kiện biên Neumann đơn giản là ∂u(x, 0) / ∂y = ∂u(x, 1) / ∂y = 0 được áp dụng tại y = 0 và y = 1. Dĩ nhiên là bài toán này có thể giải được theo cách giải tích để cho ta

u(x, y) = (1 − x2) (2 y3 − 3 y2 + 1). 

Các Hình [p2dc] và [p2dd] cho thấy sự so sánh giữa các nghiệm giải tích và sai phân hữu hạn cho N = J = 64. Có thể thấy rằng nghiệm sai phân hữu hạn khá trùng khớp với nghiệm giải tích.

p2dc

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trụcx tại y = 0. 5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2dc]

p2dd

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trụcy tại x = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2dd]

Ví dụ tính toán trường tĩnh điện 2 chiều


Ta hãy làm một ví dụ tính toán trường tĩnh điện 2 chiều. Xét một dây dẫn mang điện chạy song song trục của một đường kênh dẫn rỗng hình hộp chữ nhật. Giả sử rằng các đỉnh của hình hộp này nằm ở (x, y) = (0, 0), (0, 1), (1, 0), và (1, 1). Cũng giả sử rằng dây dẫn mang điện tích đều bằng 1 đơn vị trên mỗi đơn vị dài của dây. Điện thế φ(x, y) bên trong kênh dẫn thỏa mãn [xem PT ([estat])]

∂2φ(x, y) / ∂x2 + ∂2φ(x, y) / ∂y2 = v(x, y) =  − δ(x − x0) δ(y − y0),   (estat1)

troing đó (x0, y0) là các tọa độ của dây. Ở đây để cho tiện ta đã chuẩn hóa các đơn vị sao cho hệ số ε0 được hấp thụ trong quá trình chuẩn hóa. Giả sử rằng hộp được nối đất, điện thế sẽ tuân theo các điều kiện biên Dirichlet φ = 0 tại x = 0, x = 1, y = 0, và y = 1. Ta cần tìm nghiệm trong vùng 0 ≤ x ≤ 1 và0 ≤ y ≤ 1.

coulomb1

Hình. Đồ thị đường đồng mức của điện thế hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại(x, y) = (0. 5, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub1]

Lưu ý rằng khi làm rời rạc PT ([estat1]) vế phải trở thành

v(xi, yj) =  − 1 / (δx δy)

tại điểm nút lưới sát dây dẫn nhất, với v(xi, yj) = 0 trên các điểm nút lưới còn lại. Ở đây, δx và δy lần lượt là các khoảng cách nút lưới theo các phương xvà y.

coulomb2

Hình. Đồ thị véc-tơ thể hiện hướng của điện trường hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 5, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub2]

Các Hình [coub1] và [coub2] biểu diễn điện thế φ(x, y) và điện trường E =  − ∇φ hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật, tức là tại (x0, y0) = (0. 5, 0. 5). Việc tính toán được thực hiện với chương trình giải bài toán Poisson 2 chiều nêu trên với N = J = 128.

coulomb3

Hình. Đồ thị đường đồng mức của điện thế hình thành do một dây dẫn được tích điện đặt lệch tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại(x, y) = (0. 25, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub3]

Các Hình [coub3] và [coub4] biểu diễn điện thế φ(x, y) và điện trường E =  − ∇φ hình thành do một dây dẫn được tích điện đặt lệch khỏi tâm kênh dẫn chữ nhật, tức là tại (x0, y0) = (0. 25, 0. 5). Việc tính toán được thực hiện với chương trình giải bài toán Poisson 2 chiều nêu trên với N = J = 128.

coulomb4

Hình. Đồ thị véc-tơ thể hiện hướng của điện trường hình thành do một dây dẫn được tích điện đặt lệch tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 25, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub4]

Bài toán 3 chiều

Những kĩ thuật được thảo luận trong các Mục [fftsine] và [fftcos] để giải phương trình Poisson trong không gian hai chiều với một lớp các điều kiện biên hạn hẹp có thể được dễ dàng mở rộng ra cho không gian 3 chiều. Trong trường hợp 3 chiều, ta cần phải thực hiện biến đổi Fourier theo hai chiều (chẳng hạn theo các phương y và z) để giản hóa bài toán về một ma trận ba đường chéo nhưng các phương trình không ràng buộc nhau. Những phương trình này có thể được nghịch dảo theo cách thông thường, và nghiệm được tái lập lại bằng biến đổi Fourier ngược hai lần.

  1. Xem Numerical recipes in C: the art of scientific computing, W.H. Press, S.A. Teukolsky, W.T. Vettering, and B.R. Flannery (Cambridge University Press, Cambridge, England, 1992), Mục 12.2.↩

  2. Xem http://www.fftw.org↩

  3. Điều này đúng cho thư viện phiên bản 2 (dùng cho khóa học này), nhưng phiên bản 3 thì khác.

Tags

Post a Comment

إرسال تعليق

أحدث أقدم