Biến đổi Hartley rời rạc và họ sin, cosin

Biến đổi Hartley rời rạc và họ sin, cosin

• Định nghĩa biến đổi Hartley tích phân với hạt nhân cas, nối nó với biến đổi Fourier bằng S=e-io và \psi=\mathrm{Re}\,S-\mathrm{Im}\,S.
• Tính DHT của các chuỗi ngắn, dùng tính đối ngược (cùng công thức cho biến đổi thuận và nghịch) và các định lý dịch, đảo, sai phân, nội dung bình phương.
• Thực hiện chập bằng DHT (hai phép nhân thực thay bốn) và lấy phổ công suất trực tiếp từ DHT.
• Cài đặt thuật toán Hartley nhanh: hoán vị đảo bit, công thức phân rã chẵn lẻ, và đối chiếu số phép nhân với FFT.
• Phân biệt DFT, DHT, DST, DCT-1 và DCT-2, và giải thích vì sao DCT-2 được dùng để nén ảnh.

Dùng nút, chấm tròn hoặc phím ← → để chuyển slide.

Một biến đổi thực, đối ngược chặt chẽ

Cho dạng sóng thực V(t), ta định nghĩa \psi(\omega)=(2\pi)^{-1/2}\int V(t)(\cos\omega t+\sin\omega t)\,dt. Rõ ràng \psi là tổng các biến đổi cosin và sin hai phía, nên nghịch đảo có cùng dạng: V(t)=(2\pi)^{-1/2}\int\psi(\omega)(\cos\omega t+\sin\omega t)\,d\omega. Các hệ thức này do Hartley đưa ra năm 1942 và không bao giờ biến mất khỏi tài liệu kỹ thuật nhưng ít được biết. Thuận và nghịch giống hệt nhau, và nếu V(t) thực thì \psi(\omega) cũng thực (Bracewell, tr. 293). Trong chương này Bracewell bỏ thừa số (2\pi)^{-1/2} và chuyển sang ký hiệu quen thuộc với 2\pi\times tần số.

Từ Hartley sang Fourier: phản xạ và cộng

Viết \psi(\omega)=e(\omega)+o(\omega) với e,o là phần chẵn và lẻ: e(\omega)=(2\pi)^{-1/2}\int V\cos\omega t\,dt và o(\omega)=(2\pi)^{-1/2}\int V\sin\omega t\,dt. Khi đó biến đổi Fourier S(\omega)=(2\pi)^{-1/2}\int V(\cos\omega t-i\sin\omega t)dt=e(\omega)-io(\omega) (3): từ \psi ta lấy S bằng phản xạ và cộng đơn giản. Ngược lại \psi(\omega)=\mathrm{Re}\,[S(\omega)]-\mathrm{Im}\,[S(\omega)] (4): từ S ta được \psi là tổng phần thực và phần ảo đổi dấu (tr. 293 đến 294).

Ví dụ: V(t)=e^{-t} với t>0

Bỏ thừa số (2\pi)^{-1/2}: biến đổi Fourier S(\omega)=\dfrac{1-i\omega}{1+\omega^2} và biến đổi Hartley \psi(\omega)=\dfrac{1+\omega}{1+\omega^2} (Bracewell, tr. 294, hình 12.1); \psi thực và rõ ràng không đối xứng. Tích phân số \int_0^\infty e^{-t}\mathrm{cas}\,\omega t\,dt cho \psi(0) = 1, \psi(1) = 1, \psi(2) = 0.6 và \psi(-1) = 0. Tại \omega=1: \mathrm{Re}\,S = 0.5 và \mathrm{Im}\,S = -0.5, nên \psi=\mathrm{Re}-\mathrm{Im} = 1. Từ phần chẵn và lẻ của \psi ta khôi phục S, đúng như phép dựng trong sách.

Tính đối ngược bằng số

Với V(t)=e^{-\pi(t-0.5)^2} (Gauss lệch tâm, không chẵn không lẻ) và hạt nhân \mathrm{cas}\,2\pi st: biến đổi Fourier là e^{-\pi s^2}e^{-i\pi s}, nên \psi(s)=e^{-\pi s^2}[\cos\pi s+\sin\pi s]. Áp dụng chính công thức đó lần thứ hai cho \psi trả lại V: tại t=0.3 ta có V = 0.8819 và giá trị tái tạo lệch 1e-16. Không cần dấu i hay cặp \pm như Fourier; đó là điều Hartley gọi là đối xứng hơn.

Cái nhìn trên mặt phẳng phức

Vẽ F(s)=R+iI như quỹ tích trên mặt phẳng phức (hình 2.10). Xoay hệ trục đi -45^\circ: R'(s)+iI'(s)=[R(s)+iI(s)]e^{i\pi/4}. Thành phần theo hướng R' chính là biến đổi Hartley (nhân 1/\sqrt2), còn thành phần S' không thêm thông tin nếu f thực. Quỹ tích khôi phục được từ R và I cùng nhau; sự vắng dư thừa do đối xứng của R và I cho phép khôi phục f thực chỉ từ một hướng R', hay từ mọi hướng trừ hai trục R và S (tr. 307). Số đo tại \omega=1 của ví dụ: R' = 0.7071, và \psi/R' = 1.4142. Với dữ liệu phức, hướng S' mang thông tin còn lại và biến đổi Hartley trở nên phức (Millane, 1994).

Ý nghĩa vật lý: pha từ biên độ

Trường bức xạ của ăng ten là biến đổi Fourier của phân bố mặt mở, và thấu kính đổi phân bố này thành phân bố khác biểu diễn được tương tự, nên có vẻ Hartley chỉ là cấu trúc toán học. Nhưng có quan hệ mật thiết: mọi tính chất của biến đổi này đều có tính chất tương ứng ở biến đổi kia. Thiết bị đầu tiên dùng bộ tách chùm quang và giao thoa kế Michelson (Bracewell và cộng sự, 1985); sau đó thay bằng phần tử vi ba để mở rộng bước sóng (Villasenor và Bracewell, 1988). Trong mặt phẳng Hartley pha được mã hóa chỉ bằng biên độ; điều đáng quan tâm vì ở tử ngoại và tia X pha quan trọng nhưng khó đo (tr. 307).

Ralph Hartley và lịch sử

Hartley là học giả Rhodes, thành viên Viện Kỹ sư Vô tuyến, phụ trách nghiên cứu đường dây điện thoại tại Bell Telephone Laboratories (1918 đến 1929); nhiều thế hệ kỹ sư điện biết tên ông qua bộ dao động Hartley, ví dụ giáo khoa về nguồn sóng sin dùng một triode có hồi tiếp lấy từ cuộn cảm có đầu ra. Số đặc biệt của Proc. IEEE (3/1994, Olejniczak và Heydt) có hàng trăm bài. Khi thuật toán Hartley nhanh xuất hiện năm 1984 có nghi ngờ về ưu điểm tốc độ, thực và khả đảo; nhưng hàng trăm bài báo ra đời, chip được sản xuất, và ý nghĩa vật lý được chứng minh trong phòng thí nghiệm với vi ba và ánh sáng. Sớm rõ rằng quan hệ với Fourier là tương đương hoàn toàn (Bracewell, tr. 296 đến 297, 308 đến 309).

Tự kiểm tra phần 1

• Biến đổi Hartley của e^{-t}H(t) tại \omega=-1 bằng bao nhiêu?
• Làm thế nào lấy biến đổi Fourier từ biến đổi Hartley?
• Vì sao biến đổi Hartley của hàm thực là thực và khả đảo?

Gợi ý: (1+\omega)/(1+\omega^2)=0; lấy phần chẵn trừ i nhân phần lẻ; vì hạt nhân cas thực và nghịch đảo cùng dạng.

Định nghĩa DHT

Với biến rời rạc \tau=0,\ldots,N-1 và f(\tau) thực, H(\nu)=N^{-1}\sum_{\tau=0}^{N-1}f(\tau)\,\mathrm{cas}(2\pi\nu\tau/N), trong đó \mathrm{cas}\,\theta=\cos\theta+\sin\theta. So sánh với DFT F(\nu)=N^{-1}\sum f(\tau)e^{-i2\pi\nu\tau/N}. Nghịch đảo: f(\tau)=\sum_{\nu=0}^{N-1}H(\nu)\,\mathrm{cas}(2\pi\nu\tau/N) (Bracewell, tr. 295). Chứng minh dùng tính trực giao \sum_{\nu}\mathrm{cas}(2\pi\nu\tau/N)\,\mathrm{cas}(2\pi\nu\tau'/N)=N nếu \tau=\tau' và 0 nếu khác. Số đo với N=16: độ lệch của tổng trực giao là 5e-14, và biến đổi rồi nghịch đảo dữ liệu ngẫu nhiên trả lại dữ liệu với độ lệch 6e-15.

Bán đối xứng và bản đối xứng chặt

DHT là bán đối xứng (quasi-symmetrical): thuận và nghịch chỉ khác thừa số N^{-1} quen thuộc từ DFT, và nó thực. Để đối xứng chặt có thể đổi N^{-1} thành N^{-1/2} ở cả hai chiều, như Hartley làm với tích phân. Nhưng khi tính toán thường tiện dừng ở việc nhân hằng số, để dồn vào thừa số chuẩn hóa hay tỉ lệ đồ họa sau này. Với N^{-1/2}, ma trận \mathbf C_{\nu\tau}=N^{-1/2}\mathrm{cas}(2\pi\nu\tau/N) đối xứng và \mathbf C^2=\mathbf I: số đo N=8, \|\mathbf C^2-\mathbf I\| = 2e-15. Nghĩa là \mathbf C là ma trận đối hợp (involution) và mọi trị riêng của nó là \pm1 (tr. 296).

Từ DHT sang DFT và ngược lại

Tách H(\nu)=E(\nu)+O(\nu), với phần chẵn E(\nu)=\tfrac12[H(\nu)+H(N-\nu)] và phần lẻ O(\nu)=\tfrac12[H(\nu)-H(N-\nu)]. Khi đó F(\nu)=E(\nu)-iO(\nu); ngược lại H(\nu)=\mathrm{Re}\,F(\nu)-\mathrm{Im}\,F(\nu) (tr. 296). Phổ công suất và các quan hệ Hermite theo đó đọc thẳng được từ H. Số đo với dữ liệu ngẫu nhiên N=16: F=E-iO lệch khỏi FFT 8e-16, và H=\mathrm{Re}\,F-\mathrm{Im}\,F lệch 1e-15. Muốn có DHT nhanh nhất từ phần mềm sẵn có: chạy FFT phức rồi lấy phần thực trừ phần ảo, chỉ dùng nửa các giá trị (tr. 324).

Các cặp ngắn: N=2, N=4, N=8

N=2: \{a\ b\} có DHT \tfrac12\{a+b\ \ a-b\}; với \{1,2\} ra 1.5 và -0.5 (tr. 309). N=4: \{1,2,3,4\} có DHT \{2.5,-1,-0.5,0\} (tr. 310). N=8: \{1,2,3,4,5,6,7,8\} có DHT \{4.5,\,{-1.7071},\,{-1},\,{-0.7071},\,{-0.5},\,{-0.2929},\,0,\,0.7071\} (tr. 314; sách in làm tròn): các số đo là H(1) = -1.7071, H(3) = -0.7071, H(5) = -0.2929, H(7) = 0.7071, và H(6) = 0. Kiểm tra: tổng dữ liệu bằng N\times H(0) = 36, và tổng các giá trị DHT bằng f(0) = 1.

Ví dụ 16 điểm: mũ cụt

Lấy f(0)=0.5 (trung bình hai giá trị tại điểm gián đoạn của V(t)) và f(\tau)=e^{-\tau/2} với 1\le\tau\le15: biểu diễn hàm liên tục e^{-t/2}H(t) bằng N=16 mẫu (hình 12.2). DHT rất giống các mẫu của biến đổi liên tục lấy cách 1/16; khác biệt nhỏ do cắt cụt và chồng phổ, y như với DFT. Số đo: H(0) = 0.1275, H(1) = 0.1386 so với giá trị liên tục \psi(1/16)/16 = 0.138, với \psi(s)=(0.5+2\pi s)/(0.25+4\pi^2s^2). Đó là bản Hartley của bài toán DFT so với FT ở module 14.

Ví dụ nhị thức: cái gù trơn

Dãy nhị thức 1,6,15,20,15,6,1 biểu diễn mẫu của xung trơn; để được kết quả đơn giản nhất đặt đỉnh tại \tau=0: f=\{20,15,6,1,0,\ldots,0,1,6,15\} (N=16, chẵn quanh 0). Do đó DHT bằng DFT thực: H(\nu)=\{4,\ 3.56,\ 2.49,\ 1.32,\ 0.5,\ 0.12,\ 0.01,\ 0,\ \ldots\} và đối xứng H(N-\nu)=H(\nu): xung trơn đỉnh tại \nu=0 (hình 12.3). Các số đo: H(1) = 3.56, H(2) = 2.49, H(4) = 0.5. (Bản sao của sách dùng ở đây in số 0 ở vị trí \nu=4, trong khi tính đối xứng H(12)=0.5 đòi hỏi 0.5; có thể là lỗi quét hoặc in.)

Vì sao N số thực thay được 2N số thực

Thoạt nhìn lạ: DFT có N số phức, tổng 2N số thực. Nhưng tính chất Hermite làm dư một nửa: N/2 số thực đủ cho biến đổi cosin và N/2 cho biến đổi sin, hợp thành N hệ số DHT không suy biến do đối xứng. Hàm cas là sóng sin dịch 45^\circ, nên đáp ứng đồng đều với thành phần cosin và sin. Nếu dùng hạt nhân \sqrt2\sin(\theta+\alpha) với \alpha tùy ý thì đáp ứng không đều nhưng không mất thông tin trừ khi \alpha=0,\pi/2,\ldots; nghịch đảo dùng một tổ hợp của \sin\theta và \cos\theta với các trọng số phụ thuộc \alpha (Bracewell, tr. 298; các trọng số 1/(\sqrt2\sin\alpha) và 1/(\sqrt2\cos\alpha) dưới đây là kết quả của chúng tôi). Với \alpha=0.6: khôi phục dữ liệu lệch 5e-15.

Tự kiểm tra phần 2

• DHT của \{1,2,3,4\} là gì?
• Công thức nghịch đảo của DHT khác công thức thuận ở đâu?
• Vì sao DFT của chuỗi thực có dư một nửa thông tin?

Gợi ý: \{2.5,-1,-0.5,0\}; chỉ ở thừa số N^{-1}; vì F(N-\nu)=\overline{F(\nu)}.

Các định lý khớp hẳn

Có một định lý Hartley cho mỗi định lý Fourier; một số khớp đúng: \sum H(\nu)=f(0) và \sum f(\tau)=NH(0). Đảo: f(-\tau) có DHT H(-\nu), với chỉ số ngoài 0\ldots N-1 hiểu modulo N. Nội dung bình phương (Parseval): \sum|f(\tau)|^2=N\sum|H(\nu)|^2, giống DFT nhưng không có liên hợp phức, chỉ có số thực (Bracewell, tr. 300 đến 301). Số đo với \{1,\ldots,8\}: \sum f^2 = 204 và N\sum H^2 = 204; độ lệch của định lý đảo với dữ liệu ngẫu nhiên là 1e-15. Kiểm tra bằng DFT: \sum|F|^2=\tfrac12\sum[H(\nu)^2+H(N-\nu)^2] trùng với \sum H^2.

Định lý dịch và sai phân

Khác Fourier, dịch là: DHT của f(\tau+a) bằng H(\nu)\cos(2\pi\nu a/N)-H(-\nu)\sin(2\pi\nu a/N) (tr. 300), cần cho thuật toán Hartley nhanh. Định lý sai phân bậc nhất: f(\tau+1)-f(\tau) có DHT [\cos(2\pi\nu/N)-1]H(\nu)-\sin(2\pi\nu/N)H(N-\nu) (tr. 301). Số đo với N=16 ngẫu nhiên, a=3: định lý dịch lệch 2e-15, và định lý sai phân lệch 1e-15. Chú ý dấu: hoán đổi dấu trừ thành cộng cho sai số cỡ 0.7, nên định lý phải viết đúng dấu.

Định lý chập tổng quát

Nếu f=f_1*f_2 là chập vòng thì H(\nu)=H_1(\nu)H_{2e}(\nu)+H_1(-\nu)H_{2o}(\nu), trong đó H_2=H_{2e}+H_{2o} là phần chẵn và lẻ (tr. 299; như định lý chập của DFT, f_1*f_2\supset NF_1F_2, còn có thừa số N khi H chứa N^{-1}). Số đo với hai dãy ngẫu nhiên N=16: hệ thức đúng với thừa số N, độ lệch 1e-14. Nếu f_2 chẵn thì H_{2o}=0 và H=NH_1H_2: chỉ một tích thực mỗi phần tử. Chập bằng DFT cần một phép nhân phức mỗi phần tử tức 4 phép nhân thực; theo DHT chỉ 2 trong trường hợp tổng quát.

Quy trình chập bằng DHT

Trong xử lý ảnh và lọc số thường một hàm là chẵn (bộ lọc đối xứng). Khi đó: f_1*f_2=\mathrm{DHT}[(\mathrm{DHT}\,f_1)\times(\mathrm{DHT}\,f_2)], tức lấy hai DHT, nhân từng số hạng hai dãy thực, và lấy một DHT nữa (thêm thừa số N tùy quy ước chuẩn hóa). Trước đó kéo dài hai dãy gấp đôi bằng số 0 để có chỗ cho kết quả; số 0 có thể ở đầu, cuối hay hai bên. Nếu một dãy ngắn hơn nhiều, như lọc số thông thường, thì chập trực tiếp được khuyên dùng. Với \{1,2,3,4\} và bộ lọc chẵn \{0.25,0.5,0.25\} độn đến 16: kết quả khớp chập trực tiếp với độ lệch 4e-15; chuỗi thu được \{0.25,\ 1,\ 2,\ 3,\ 3.25,\ 1\} có phần tử lớn nhất 3.

DHT hai chiều

Ảnh f(\tau_1,\tau_2) dạng ma trận M\times N có DHT là ma trận thực H(\nu_1,\nu_2)=M^{-1}N^{-1}\sum\sum f\,\mathrm{cas}[2\pi(\nu_1\tau_1/M+\nu_2\tau_2/N)] và nghịch đảo cùng công thức không có thừa số (tr. 299). Nó tách thành N biến đổi một chiều. Lợi ích: có thể trình bày trọn vẹn phân tích phổ ảnh bằng một sơ đồ (Fourier cần hai: phần thực và ảo, hoặc biên độ và pha); lưu trữ một mảng thực. Phổ công suất hai chiều giống nhau ở cả hai biến đổi. Số đo với ảnh 6\times5: H=\mathrm{Re}\,F-\mathrm{Im}\,F với F là DFT hai chiều, độ lệch 6e-16; và biến đổi hai lần trả lại ảnh với độ lệch 4e-15. Ba chiều: Hao và Bracewell (1987); bốn chiều: Buneman (1987).

Biến đổi cas-cas

Trong hai chiều có thể dùng hạt nhân \mathrm{cas}\,2\pi\nu_1\tau_1\,\mathrm{cas}\,2\pi\nu_2\tau_2, tách được (Perkins, 1987). Biến đổi cas-cas cũng đối ngược (Bracewell, 1983) và liên hệ với DHT chuẩn giống như chọn giữa các thành phần Fourier dạng \cos2\pi(\nu_1\tau_1+\nu_2\tau_2), tức nếp gấp có đường không song song, và \cos\cos có đường không vuông góc nhau. Tính tách được không chứng tỏ là cơ sở để ưu tiên, nên hạt nhân \mathrm{cas}[2\pi(\nu_1\tau_1+\nu_2\tau_2)] thành chuẩn (tr. 300). Số đo trên ảnh 8\times8: biến đổi cas-cas hai lần trả lại ảnh (độ lệch 6e-15) nhưng khác DHT chuẩn: chênh lệch lớn nhất 1.52.

Phổ công suất trực tiếp từ DHT

Nếu muốn phổ công suất |F(\nu)|^2=R^2+X^2 thì không cần đi qua số phức: |F(\nu)|^2=\tfrac12[H(\nu)^2+H(N-\nu)^2] (tr. 314). Với \{8,7,6,5,4,3,2,1\} (Bài 4): H=\{4.5,\ 1.7071,\ 1,\ 0.7071,\ 0.5,\ 0.2929,\ 0,\ -0.7071\}, phần chẵn E(1) = 0.5, phần lẻ O(1) = 1.2071, nên F(1)=E-iO = 0.5 - 1.2071i và |F(1)|^2 = 1.7071, khớp công thức trực tiếp (2e-15). Độ lệch trên mọi \nu với dữ liệu ngẫu nhiên: 3e-16.

Bài tập của chương

Bài 2: dãy \{1,2,\ldots,2^P\} luôn có đúng một phần tử DHT bằng 0. Từ P=2 (chỉ số 3) và P=3 (chỉ số 6), quy tắc là \nu=3N/4; với N=16 chỉ số 12 và N=64 chỉ số 48. Chứng minh: H(3N/4)\propto\sum(\tau+1)[\cos(3\pi\tau/2)+\sin(3\pi\tau/2)]; nhóm bốn số hạng liên tiếp cho tổng 0. Bài 8: \{1234\} lặp bốn lần có DHT khác 0 chỉ tại \nu=0,4,8 với giá trị 2.5, -1, -0.5. Bài 10: \{-2,2,-2,2,\ldots\} là sóng Nyquist, DHT chỉ có H(4) = -2, phổ công suất tập trung tại \nu=N/2.

Tự kiểm tra phần 3

• Định lý dịch của DHT khác của DFT thế nào?
• Khi nào chập bằng DHT chỉ cần một tích thực mỗi phần tử?
• Lấy phổ công suất từ DHT như thế nào?

Gợi ý: có thêm số hạng H(-\nu)\sin; khi một hàm chẵn; \tfrac12[H(\nu)^2+H(N-\nu)^2].

Nhu cầu và các ưu điểm

Chuỗi N số thực có DHT là chuỗi cùng độ dài và cũng thực; áp dụng cùng công thức lần hai trở về dãy ban đầu. Không cần quản lý phần thực và ảo ở hai mảng, hay đan xen trong một mảng gấp đôi, và không cần nghịch đảo khác với thuận. Nếu cần phần thực và ảo của DFT thì lấy được từ phần chẵn và lẻ của DHT; phổ công suất lấy trực tiếp. Chương trình FFT một chiều cho dữ liệu thực bỏ tính các hệ số liên hợp, nhưng đầu ra phức, nên phải lưu thêm chương trình nghịch đảo riêng, và không chương trình nào nhận được đầu ra của chính nó. Bracewell nêu FHT nhanh hơn FFT phức theo một hệ số tùy hoàn cảnh (tr. 308).

Công thức phân rã chẵn lẻ

Với N=2: \{a\ b\} có DHT \tfrac12\{a+b\ \ a-b\}. Chia f thành dãy chỉ số chẵn (DHT H_1) và lẻ (DHT H_2), mỗi dãy dài N/2. Dùng định lý đồng dạng (kéo giãn: DHT lặp lại) cho phần chẵn và định lý dịch với a=-1 cho phần lẻ: H(\nu)=\tfrac12\left[H_1(\nu)+H_2(\nu)\cos\tfrac{2\pi\nu}N+H_2(N-\nu)\sin\tfrac{2\pi\nu}N\right] với H_1,H_2 hiểu tuần hoàn chu kỳ N/2; ở tài liệu các thừa số \tfrac12 được dồn lại đến cuối (tr. 309 đến 310). Số đo: với N=16 ngẫu nhiên công thức khớp DHT trực tiếp, độ lệch 4e-15. Ở các giai đoạn đầu sau hoán vị, \cos,\sin chỉ nhận 0, 1, -1 nên chỉ còn cộng, trừ.

Hoán vị đảo bit

Hoán vị sắp xếp lại dãy để chia đôi liên tục tới các cặp: với 8 phần tử thứ tự là f(0),f(4),f(2),f(6),f(1),f(5),f(3),f(7). Viết chỉ số i nhị phân (000 đến 111), đảo các chữ số rồi đổi lại thập phân ra chỉ số mới j. Một nửa chỉ số (0, 2, 5, 7) không đổi vì biểu diễn nhị phân đối xứng, nên bỏ qua chúng tiết kiệm thời gian; hoán vị chỉ còn hai cặp đổi chỗ 1\leftrightarrow4 và 3\leftrightarrow6. Đồ thị tán xạ j theo i gồm những ô lặp lại (chú diều N=4, con ếch N=8) (tr. 321 đến 323). Số đo: số chỉ số cố định là 4 với N=8 và 32 với N=1024 (2^{\lceil P/2\rceil}); với N=8 có 2 cặp đổi chỗ.

Ví dụ tay: \{1,2,\ldots,8\}

Bảng 12.3: bước 1 tách \{1357\} và \{2468\}, bước 2 tách thành \{15\}, \{37\}, \{26\}, \{48\}; dãy hoán vị f_0=\{1,5,3,7,2,6,4,8\}. Mỗi cặp \{a\ b\} cho \{a+b\ \ a-b\} (bỏ thừa số \tfrac12): f_1=\{6,-4,10,-4,8,-4,12,-4\}, với 6=1+5, -4=1-5, 10=3+7, -4=3-7... Sau hai giai đoạn ghép với \cos,\sin của tám phần của vòng, ra 8H=\{36,\ -13.657,\ -8,\ -5.657,\ -4,\ -2.343,\ 0,\ 5.657\}; chia 8 cho H=\{4.5,-1.7,-1,-0.7,-0.5,-0.3,0,0.7\} (tr. 311 đến 314). Số đo: f_1(0) = 6, f_1(1) = -4, 8H(1) = -13.6569. Kiểm tra: tổng dữ liệu bằng 8 lần phần tử đầu của DHT.

Đếm phép nhân: FHT so với FFT

Cài đặt đệ quy trực tiếp công thức phân rã và đếm các phép nhân với hệ số khác 0,\pm1. Với N=1024 số phép nhân thực của FHT là 14344; trong FFT chia đôi có 3586 phép nhân phức thực sự, tức 14344 phép nhân thực (bốn thực cho một phức), bằng nhau ở dạng thô này: lợi thế của FHT nằm ở chỗ nó tính N số thực (nửa bộ nhớ và đầu ra) thay vì N số phức, và dữ liệu vào thực không cần dùng nửa phổ Hermite dư. Bracewell báo cáo bằng thực nghiệm trên máy tính cá nhân: hệ số Y nằm giữa 3.1 và 5.8 đối với chương trình FFT dịch, khoảng 1.3 nếu dùng FHT để lấy Fourier (tr. 315): số của thời gian chạy trên một máy, không tái lập ở đây.

Giản đồ vằn: thời gian theo thành phần

Bracewell phân tích thời gian chạy theo N trong sơ đồ vằn (hình 12.5). Ba loại thành phần: chi phí cố định (overhead) chết dần khi chuẩn hóa với NP; các phép tính trước sin và cos, tỉ lệ với N, nên vằn của chúng mỏng dần như P^{-1}; và hoán vị, N phép gán lặp P lần, tỉ lệ NP, nên vằn không chết dần: mọi nỗ lực tinh gọn hoán vị đều đáng giá (như thay J/2 bằng 0.5\times J). Vùng dày không đổi có độ dày 0.038 s trên máy của tác giả, nên thời gian giới hạn là 0.038NP; ở P\approx10 khoảng cao hơn một phần ba. Các giai đoạn 1 và 2 dùng phương trình ngắn của bảng 12.2 (không lượng giác) để nhanh (tr. 315 đến 317). Thuật toán cơ số 4 tiết kiệm khoảng 25 phần trăm thời gian, cho N là lũy thừa của 4 (16, 64, 256, ...) (tr. 323).

Biểu diễn ma trận và phân tích DFT

Có thể viết \mathbf H=\mathbf L_P\cdots\mathbf L_2\mathbf L_1\,\mathbf P\,\mathbf f (với \mathbf P là ma trận hoán vị và \mathbf L_i toán tử ma trận thưa của mức i), và bước sang DFT \mathbf F=\mathbf W\mathbf f bằng nhân với ma trận (phức) chuyển; kết hợp cho một cách phân tích mới của ma trận DFT \mathbf W với W=e^{-i2\pi/N} (tr. 317 đến 320). \mathbf L_1 và \mathbf L_2 chỉ có 2N phần tử khác 0 so với 3N ở các mức cao hơn; \mathbf L_1 có phần tử sin nằm trên đường chéo dốc 45^\circ. Số đo: trong ma trận cas 16\times16 có 224 phần tử khác 0 trên tổng 256 (các phần tử bằng 0 đúng là nơi \mathrm{cas}=0).

Tự kiểm tra phần 4

• Thứ tự hoán vị đảo bit của 8 phần tử là gì?
• Bao nhiêu chỉ số không đổi khi hoán vị 8 phần tử?
• Vì sao hoán vị là thành phần không chết dần trong sơ đồ vằn?

Gợi ý: 0,4,2,6,1,5,3,7; 4 (chỉ số 0, 2, 5, 7); vì tỉ lệ NP.

Biến đổi cosin dạng DFT không khả đảo

Biến đổi cosin rời rạc \sum f(\tau)\cos(2\pi\nu\tau/N) (và tương tự sin) không nghịch đảo được để khôi phục f(\tau) vì các thành phần \sin(2\pi\nu\tau/N) bị tích phân mất không để lại dấu. Chúng khả đảo khi giới hạn ở hàm f chẵn (f(\tau)=f(N-\tau) với \tau=1\ldots N-1), và biến đổi cosin/sin tích phân khả đảo nếu f(x) bằng 0 với x âm (Bracewell, tr. 301 đến 302). Số đo: chuỗi \{1,2,3,4,5,6,7,8\} và bản đảo của nó cho cùng biến đổi cosin (chênh 2e-14) dù là hai dãy khác nhau (chênh dữ liệu 6).

DST và DCT-1 khả đảo với tần số nửa nguyên

Đổi tần số thành 0,\tfrac1{2N},\tfrac1N,\tfrac3{2N},\ldots: tần số cơ bản là 1/2N, chu kỳ cơ bản 2N. DST: F_s(\nu)=\tfrac2N\sum_{\tau=1}^{N-1}f(\tau)\sin(\pi\nu\tau/N), nghịch đảo f(\tau)=\sum_{\nu=1}^{N-1}F_s(\nu)\sin(\pi\nu\tau/N); nếu áp dụng hai lần tổng \sum f\sin ta được \tfrac12Nf(\tau). Giá trị f(0) không để lại dấu vì \sin bằng 0 ở đó. DCT-1 (Rao và Yip, 1990) đối xứng chặt với N+1 giá trị: F(\nu)=\sqrt{2/N}\,k_\nu\sum_{\tau=0}^Nk_\tau f(\tau)\cos(\pi\nu\tau/N), k=1/\sqrt2 ở hai đầu và 1 ở nơi khác (tr. 301 đến 302). Số đo N=16: nghịch đảo DST lệch 4e-15; áp dụng hai lần cho hệ số 8; DCT-1 áp dụng hai lần lệch 2e-15.

Bài toán biên: dây đàn

Một dây đàn dài L bị kéo lệch tại điểm không ở giữa: điều kiện biên y(0)=y(L)=0 làm các số hạng cosin không cần, nên chỉ cần \sin(\pi x/L),\sin(2\pi x/L),\ldots để tổng hợp hình dạng dây. Nếu thay bằng [dy/dx]_{x=0}=0 thì chỉ cần cosin. Làm bằng Fourier đầy đủ vẫn được nhưng không thông dụng vì cho cả hai loại thành phần (Bracewell, tr. 302 đến 303). Ví dụ của chúng tôi (dây gảy tại a=0.3, L=1): y(x)=x/a với x<a và (1-x)/(1-a) với x>a, hệ số sin b_n=\dfrac{2\sin n\pi a}{n^2\pi^2a(1-a)}: b_1 = 0.7807, b_2 = 0.2294; tổng 300 số hạng cho y(0.5) = 0.7143 so với đúng 0.7143, lệch 8e-08. Sóng trong kênh có hai đầu đóng phân tích chỉ bằng cosin.

Các phương trình vi phân và số chiều

Các lĩnh vực vật lý về phương trình sóng, phương trình Laplace, phương trình khuếch tán và các phương trình cơ bản khác ở hai, ba và bốn chiều cho phạm vi ứng dụng rộng cho phân tích số bằng biến đổi sin và cosin rời rạc (tr. 303). Lý do: các điều kiện biên tự nhiên chọn hàm cơ sở, và nếu hàm nền sin là nghiệm của bài toán riêng thì phép biến đổi làm rời rạc hóa các toán tử vi phân thành phép nhân (module 13 và 19: nghiệm riêng của \phi''+\lambda\phi=0 là \sin k\pi t). Một minh họa số: ma trận sai phân bậc hai 50\times50 có trị riêng 2-2\cos(k\pi/51) và vectơ riêng \sin(k\pi j/51): sai lệch giữa trị riêng số và công thức 2e-15.

Nén dữ liệu và mã hóa biến đổi

Ảnh triệu mẫu, mỗi mẫu 256 mức: giảm bốn lần số bit nếu thay khối 2\times2 điểm bằng một điểm trung bình, nhưng làm mờ ảnh, mô tả bằng chập rồi lấy mẫu thô hơn. Thực nghiệm: nhiều ảnh có ít nội dung ở tần số không gian cao, gợi ý truyền biến đổi Fourier của ảnh nhưng triệt hay làm nhẹ tần số cao (chẳng hạn lượng tử hóa với ít mức hơn), rồi nghịch đảo ở đầu thu: gọi là mã hóa biến đổi (transform encoding). Còn có thể được lợi thêm bằng mã hóa một chiều (run-length, TIFF, delta-row) và phân vùng ảnh theo thống kê đồng nhất; văn bản, ASCII, thoại số ít lợi (Bracewell, tr. 303 đến 304). Số đo: nhóm 2\times2 giảm số bit 4 lần; và giảm 32 lần thêm do báo in chỉ có khoảng 8 mức xám: 256/8 = 32.

DCT-2: biến đổi của nén ảnh

Ahmed, Natarajan và Rao (1974) đưa ra {}^{\mathrm{II}}F(\nu)=\sum f(\tau)\cos[\pi\nu(2\tau+1)/2N] (thang chuẩn hóa \sqrt{2/N}\,a(\nu) với a(0)=1/\sqrt2), và nó vào chuẩn quốc tế (JPEG). DCT-2 phân tích thành tần số 0,\tfrac1{2N},\tfrac1N,\ldots như DCT-1 nhưng mẫu dịch nửa đơn vị, nên hàm cơ sở không còn thuần cosin; các vectơ cơ sở vẫn trực giao. Bracewell nhấn mạnh biến đổi này không khả đảo chặt (áp dụng hai lần không trả lại vật): số đo \|\mathbf C^2-\mathbf I\| = 1.99 nhưng \mathbf C\mathbf C^T=\mathbf I lệch 2e-15. Các định lý dịch và chập trở nên mờ đục. Ưu điểm được nêu: hiệu quả thực nghiệm tốt trên ảnh có thống kê xấp xỉ Markov bậc nhất đồng đều (tr. 304).

Vì sao DCT-2 hợp Markov bậc nhất: gom năng lượng

Quá trình Markov bậc nhất có tự tương quan r(k)=\rho^{|k|} (module 14, AR(1)). Với \rho=0.95 và khối N=16, xếp các phương sai hệ số (đường chéo của \mathbf TR\mathbf T^T) giảm dần và lấy bốn hệ số lớn nhất: phần năng lượng giữ lại là 0.9557 với DCT-2, 0.9364 với DHT chuẩn hóa, 0.9559 với biến đổi Karhunen-Loève (giới hạn tối ưu, module 13), và chỉ 0.25 nếu không biến đổi (lấy 4 mẫu). Vậy DCT-2 gần như đạt giới hạn tối ưu; đó là lý do thực nghiệm nó được ưa chuộng để nén (kiểm chứng số của nhận xét ở tr. 304; giá trị chính xác là kết quả của chúng tôi).

Ví dụ tần số âm và chỉ số DHT

Bracewell nhắc rằng tần số f nằm ở hai chỗ \nu=\pm Nf mod N; thành phần cosin có giá trị bằng nhau, thành phần sin bằng nhau nhưng trái dấu (tr. 305 đến 306). Với N=8: \{1,0,-1,0,1,0,-1,0\} (\cos2\pi\tau/4) có DHT \tfrac12 tại \nu=2 và 6: 0.5 và 0.5; \{0,1,0,-1,0,1,0,-1\} (\sin) có DHT 0.5 tại \nu=2 và -0.5 tại \nu=6; thêm dc 2: \{3,2,1,2,3,2,1,2\} có DHT 2 tại \nu=0. Với các biến đổi khác, DCT-1 và DCT-2 của chính những chuỗi này cho số khác, vì chúng không phải Fourier: nên đọc phổ DCT cần cẩn thận về ý nghĩa tần số (tr. 306).

Tự kiểm tra phần 5

• Vì sao biến đổi cosin dạng DFT không khả đảo?
• DST khả đảo bằng cách nào?
• Vì sao DCT-2 phù hợp dữ liệu Markov bậc nhất?

Gợi ý: thành phần sin bị tích phân mất; dùng tần số 0,\tfrac1{2N},\tfrac1N,\ldots và bỏ \tau=0; gom gần hết năng lượng vào vài hệ số như KLT.

Điều cần mang theo

• Hartley: \psi=\mathrm{Re}\,S-\mathrm{Im}\,S với hạt nhân \mathrm{cas}=\cos+\sin; thuận và nghịch cùng dạng; là thành phần của quỹ tích Fourier theo trục quay 45^\circ.
• DHT: N số thực vào, N số thực ra; F=E-iO, H=\mathrm{Re}F-\mathrm{Im}F; \{1,2,3,4\}\to\{2.5,-1,-0.5,0\}.
• Định lý: dịch H\cos-H(-\nu)\sin; Parseval \sum f^2=N\sum H^2; chập H_1H_{2e}+H_1(-\nu)H_{2o}; phổ công suất \tfrac12[H(\nu)^2+H(N-\nu)^2].
• FHT: hoán vị đảo bit rồi phân rã chẵn lẻ H=\tfrac12[H_1+H_2\cos+H_2(N-\nu)\sin]; số phép nhân thô bằng FFT nhưng chỉ dùng số thực.
• DST, DCT-1 khả đảo (tần số nửa nguyên); DCT-2 (JPEG) trực giao nhưng không đối hợp, gom năng lượng của quá trình Markov bậc nhất gần như Karhunen-Loève.