Tích chập, tích chập nối tiếp và tương quan

Tích chập, tích chập nối tiếp và tương quan

• Định nghĩa tích chập và hình dung nó bằng bốn cách: diện tích, làm mượt, chồng chất, gập.
• Tính tích chập của hàm mũ cắt cụt và của các dãy số bằng tay, kiểm bằng tổng các số hạng.
• Đảo tích chập nối tiếp (chia dãy), dùng bảng dãy nghịch đảo, và viết tích chập dưới dạng ma trận.
• Tính tích chập bằng máy, biết khi nào FFT bị cuộn vòng.
• Dùng tự tương quan, tương quan chéo, tương quan bậc ba và phổ năng lượng, và biết chúng bỏ mất thông tin gì.

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

Một khái niệm nhiều tên

Tích chập xuất hiện ở nhiều nơi với nhiều tên: Faltung (tiếng Đức), tích thành phần ("composition product"), tích phân chồng chất, tích phân Duhamel, định lý Borel, trung bình trượt (có trọng số), hàm tương quan chéo, làm mượt, làm nhòe, quét (Bracewell, tr. 24).

Định nghĩa

Tích chập của f(x) và g(x) cũng là một hàm của x, gọi là h(x) (tr. 24).

• ubiến giả để lấy tích phân
• g(x-u)g đã được đảo chiều rồi dịch tới x

Hàm của hàm: cần biết f trên cả một khoảng

Theo Volterra, h là một phiếm hàm của f. Để tính h(x_1) ta cần biết f trên cả một khoảng x, trong khi một hàm của hàm f tại x=x_1 chỉ cần f(x_1) (tr. 25). Đây là lý do tích chập mô tả "trí nhớ" của hệ.

Cách nhìn 1: diện tích của tích f(u)g(x-u)

Trong hình 3.1, tích f(u)g(x-u) được tô bóng và tung độ h(x) bằng diện tích phần tô (tr. 25). Ví dụ số: f=\Pi và g=e^{-u}H(u); tại x=1 diện tích chồng lấn là e^{-1/2}-e^{-3/2} = 0.3834 (tích phân số và công thức đóng trùng nhau).

Cách nhìn 2: làm mượt

So với f, hàm h=f*g mượt hơn, rộng hơn và có tổng biến thiên nhỏ hơn (hình 3.2, tr. 25). Ví dụ: chuỗi \pm1 đảo dấu liên tục có tổng biến thiên 48; sau khi làm trung bình trượt 5 điểm còn 11.6.

Cách nhìn 3 và 4: chồng chất và gập

Cách 3: chia f thành các cột nhỏ, mỗi cột "tan" thành một đống có dạng g đặt tại vị trí cột, rồi h(x) là tổng các đóng góp tại x (hình 3.3). Cách 4: gập g(u) quanh đường u=x/2, diện tích dưới tích f(u)g(x-u) chính là h(x) (hình 3.4) (tr. 25 đến 26).

Đại số: giao hoán, kết hợp, phân phối

Vì g bị đảo chiều trước khi nhân, tích chập giao hoán. Nó còn kết hợp (khi các tích phân tồn tại) và phân phối trên phép cộng, nên dấu * cư xử như dấu nhân trong biến đổi hình thức (tr. 26 đến 27). Notebook kiểm cả ba tính chất trên dãy ngẫu nhiên: 3 trên 3 đều đúng.

Ý nghĩa vật lý: mọi phép quan sát đều là một tích chập

Tích chập mô tả tác động của một dụng cụ khi lấy trung bình có trọng số một đại lượng vật lý trên một khoảng hẹp của một biến. Khi dạng của hàm trọng số không đổi theo vị trí, đại lượng quan sát được là tích chập chứ không phải đại lượng cần đo. Mọi phép quan sát đều bị giới hạn bởi độ phân giải của dụng cụ, chỉ riêng điều này đã làm tích chập có mặt khắp nơi (tr. 24). Sự xuất hiện của tích chập trùng với tuyến tính cộng bất biến, và với "đáp ứng điều hòa với kích thích điều hòa".

Tự kiểm tra phần 1

• Vì sao tích chập giao hoán dù g bị đảo chiều còn f thì không?
• Vì sao gọi h(x) là một phiếm hàm của f?
• Khi nào ảnh quan sát được là tích chập của vật với dụng cụ?

Gợi ý: đổi biến u\to x-u; cần f trên cả khoảng; khi hàm trọng số không đổi theo vị trí.

Tích chập của hai hàm mũ cắt cụt

Đặt E(x)=e^{-x}H(x). Với hai hằng số suy giảm dương khác nhau \alpha,\beta, Bracewell tính (tr. 27):

• \alpha,\betahằng số suy giảm

Số đo: chuỗi phân rã phóng xạ

Kết quả là hiệu của hai hàm mũ cắt cụt cùng biên độ. Nó mô tả nồng độ một đồng vị phân rã với hằng số \alpha trong khi được bổ sung từ đồng vị mẹ phân rã với hằng số \beta. Với \alpha=1,\beta=2 tại x=1: tích phân số và công thức cho 0.4651. Khi x\to\infty chỉ còn số hạng suy giảm chậm hơn (tr. 28).

Giới hạn \beta\to\alpha: bộ cộng hưởng tới hạn

Lấy giới hạn \beta-\alpha\to0 ta có \alpha E(\alpha x)*\alpha E(\alpha x)=\alpha^2xE(\alpha x). Hàm này mô tả đáp ứng của bộ cộng hưởng tắt dần tới hạn, như điện kế không dao động, với một xung tác động (tr. 28). Với \alpha=1 tại x=2: 0.2707.

Hai đuôi ngược chiều

\alpha E(-\alpha x)*\beta E(\beta x) cho một hàm nhọn tại gốc, suy giảm với hằng số \alpha về bên trái và \beta về bên phải (tr. 28 đến 29). Với \alpha=1,\beta=2: tại x=-1 là 0.2453, tại x=+1 là 0.0902, nên bên phải rơi nhanh hơn.

Dựng đồ thị: tờ giấy trượt

Để kiểm phép tính, vẽ một hàm đảo ngược lên một tờ giấy rồi trượt dọc trục hoành. Khi tờ giấy ở bên trái, tích bằng 0. Rồi tích phân bắt đầu khác 0, tăng gần tuyến tính, đạt cực đại rồi giảm. Bằng cách này ta thấy ngay các đuôi ngược chiều, chỗ không có đoạn bằng 0 và bước nhảy độ dốc tại gốc (tr. 29 đến 30). Chuyên gia làm việc này trong đầu.

Lặp tích chập xung chữ nhật

\Pi*\Pi=\Lambda (tam giác), rồi \Pi*\Pi*\Pi là một đường parabol từng khúc. Notebook lặp tích chập trên lưới mịn: tại x=0 giá trị là 0.75 (đúng 3/4) và tại x=1 là 0.125 (đúng 1/8). Mỗi lần tích chập, hình dạng mượt hơn và trông giống Gauss hơn.

Quy tắc diện tích và phương sai

Với các hàm diện tích 1, tích chập nhân diện tích (nên vẫn diện tích 1) và cộng phương sai. Ví dụ e^{-x}H(x) (phương sai 1) chập với \Pi (phương sai 1/12) cho phương sai 1.0833 và diện tích 1. Lặp n lần \Pi có phương sai n/12; độ lệch tối đa so với Gauss cùng phương sai là 0.0243 khi n=4 và 0.0092 khi n=8.

Tự kiểm tra phần 2

• Khi \beta\to\alpha kết quả tiến về hàm nào?
• Bên nào của \alpha E(-\alpha x)*\beta E(\beta x) suy giảm nhanh hơn nếu \beta>\alpha?
• Phương sai của tích chập của ba xung chữ nhật là bao nhiêu?

Gợi ý: \alpha^2xE(\alpha x); bên phải; 3/12=1/4.

Nhân đa thức chính là tích chập nối tiếp

Tích của hai đa thức a_0+a_1x+\dots và b_0+b_1x+\dots có hệ số c_0=a_0b_0, c_1=a_0b_1+a_1b_0, c_2=a_0b_2+a_1b_1+a_2b_0, ... Đó chính là các tổng của tích chập rời rạc (Bracewell, tr. 30). Notebook so polymul với convolve và thấy trùng.

Độ dài và tổng: hai phép kiểm nhanh

Tích chập nối tiếp dài hơn từng thành phần: số số hạng bằng tổng hai độ dài trừ một. Tổng các số hạng bằng tích hai tổng, một phép kiểm số rất quý. Nếu một dãy có tổng bằng 1 thì tổng tích bằng tổng dãy kia (tr. 32). Ví dụ: dãy 7 số hạng chập dãy 4 số hạng cho 10 số hạng.

Dãy nửa vô hạn và hai phía

Dãy nửa vô hạn như \{1,1\}*\{1,2,3,\dots\}=\{1,3,5,7,\dots\} vẫn cho tích chập nối tiếp hợp lệ vì mỗi số hạng chỉ có hữu hạn số hạng cộng lại. Nhưng nếu hai dãy kéo dài ngược chiều thì mỗi số hạng là tổng vô hạn và có thể không tồn tại. Với dãy hai phía phải nêu rõ gốc, ví dụ bằng mũi tên (tr. 32 đến 33).

Các dãy đặc biệt

Dãy \{\dots0,0,1,0,0,\dots\} đóng vai trò như xung: chập với nó không đổi gì. \{1,-1\} lấy sai phân bậc nhất. \frac1n\{1,\dots,1\} (n số hạng) là trung bình trượt; ví dụ trung bình tuần của số liệu ngày là \frac17\{1,1,1,1,1,1,1\}. \{1,1,1,\dots\} cho tổng chạy (tr. 33).

Thử số: tổng chạy và sai phân là hai phép nghịch đảo

Với dãy \{3,1,4,1,5,9,2,6\}, tổng chạy có số cuối 31, và lấy sai phân bằng \{1,-1\}* trả lại đúng dãy gốc. Đây là hệ thức \{1,1,1,\dots\}*\{1,-1\}=\{1,0,0,\dots\} (tr. 34): hai dãy nghịch đảo của nhau, như tích phân và đạo hàm.

Dãy nghịch đảo

Dãy nghịch đảo \{f\}^{-1} của \{f\} thỏa \{f\}^{-1}*\{f\}=\{1,0,0,\dots\}. Nếu \{f\}*\{g\}=\{h\} và biết \{f\}, \{h\} thì \{g\}=\{f\}^{-1}*\{h\}. Ví dụ \{1,2,3,4,\dots\} và \{1,-2,1\} là cặp nghịch đảo (tr. 34). Notebook kiểm: 2 trên 2 cặp đúng.

Chia dãy: tìm thừa số còn lại

Muốn giải \{1,1\}*\{g\}=\{1,3,3,1\}, có thể chia đa thức hoặc dùng cách của Bracewell: g_0=h_0/f_0, rồi mỗi số hạng tiếp theo bằng h_k trừ tổng các tích đã biết, chia cho f_0 (tr. 34 đến 36). Cả hai cho \{g\} = (1 2 1).

Bảng 3.1: các dãy và nghịch đảo

Bracewell liệt kê các cặp thường gặp (tr. 37); notebook kiểm 6 trên 6 cặp đầu bằng nhân rồi xem {1, 0, 0, ...}.

Dạng ma trận và hệ thay đổi theo thời gian

Tích chập nối tiếp viết được thành nhân một ma trận vuông (mỗi hàng là dãy \{f\} đảo ngược, dịch một bước) với một vector cột. Ma trận này là ma trận Toeplitz. Cách viết hy sinh tính giao hoán nhưng cho phép tổng quát hóa: nếu các hàng thay đổi ta được phép biến đổi tuyến tính tổng quát, không phải bất biến. Với hệ thay đổi theo thời gian (như chiết áp chạy bằng động cơ), đáp ứng không còn là tích chập và "đáp ứng điều hòa với kích thích điều hòa" cũng hỏng (tr. 37 đến 39). Notebook kiểm: nhân ma trận Toeplitz và convolve trùng nhau (yes).

Mã tính trực tiếp

Bracewell viết đoạn mã ngắn cho tích chập trực tiếp: với hai dãy độ dài l_f và l_g, tích chập có độ dài l_h=l_f+l_g-1, và mỗi h(I) cộng các tích f(J)g(K) với K=I-J+1 (tr. 40). Trong ứng dụng lọc dữ liệu dài, đoạn mã này thường chạy nhanh hơn phương pháp biến đổi viết cùng ngôn ngữ (tr. 40).

Thử số: mã tay so với convolve

Notebook dịch đoạn mã sang Python và so với numpy.convolve trên hai cặp dãy ngẫu nhiên (độ dài 200 và 37, rồi 5 và 5): kết quả trùng tới sai số làm tròn, lệch lớn nhất 5.3e-15.

Ba chế độ độ dài

Với hai dãy độ dài 10 và 4: chế độ full (mọi vị trí chồng lấn) cho 13 số hạng, same cho 10 (cùng độ dài dãy dài) và valid (chỉ nơi hai dãy chồng lấn hoàn toàn) cho 7. Khi làm mượt số liệu, hãy nhớ rằng bờ dữ liệu bị ảnh hưởng.

Đếm phép nhân: trực tiếp so với FFT

Trực tiếp cần l_fl_g phép nhân. Qua FFT cần ba phép biến đổi độ dài N\ge l_f+l_g-1, cỡ 3N\log_2N phép tính (ước lượng bậc độ lớn). Với l_f=10^5 và bộ lọc ngắn l_g=32: trực tiếp 3,200,000, FFT cỡ 6,684,672, cùng bậc. Với bộ lọc dài l_g=2000: trực tiếp 200,000,000, FFT vẫn cỡ 6,684,672, nhanh hơn rõ rệt.

Bẫy cuộn vòng

Nhân hai FFT có độ dài N nhỏ hơn l_f+l_g-1 cho tích chập vòng, không phải tích chập thẳng: phần đuôi cuộn về đầu (chương 11 sẽ gọi là tích chập tuần hoàn). Với f dài 8, g dài 3 và N=8, hai kết quả lệch nhau tới 1.09 ở các số hạng đầu. Đệm số 0 tới N\ge l_f+l_g-1 thì lệch chỉ 4.4e-16.

Ví dụ: trung bình 7 ngày loại chu kỳ tuần

Bracewell nhắc rằng trung bình trượt dùng để làm mượt số liệu khí tượng (tr. 33). Với chuỗi 364 ngày có xu hướng mùa biên độ 3 và chu kỳ tuần biên độ 2, trung bình 7 ngày (đóng vòng) đưa biên độ tuần từ 2 xuống 0 và giữ biên độ mùa 2.9982 (lý thuyết 2.9982). Hai cách tính, tích chập vòng và FFT, trùng nhau.

Thư viện ma trận và tích chập

Trong các gói phần mềm dựa trên phép toán ma trận, tích chập thường được dựng thành tích của một ma trận vuông với một vector cột, và người ta coi thời gian mất vì kém hiệu quả là không đáng kể (tr. 38). Thực hành hiện nay dùng hàm chuyên dụng (convolve, fftconvolve) vì chúng chọn cách tính hợp lý theo độ dài.

Tự kiểm tra phần 4

• Tích chập hai dãy độ dài 12 và 5 có bao nhiêu số hạng ở chế độ full và valid?
• Vì sao nhân hai FFT độ dài N=8 của dãy dài 8 và dãy dài 3 cho kết quả sai?
• Với bộ lọc rất ngắn, phương pháp nào thường kinh tế hơn?

Gợi ý: 16 và 8; vì tích chập vòng cuộn đuôi; tính trực tiếp.

Tự tương quan: định nghĩa

Tự chập của f là \int f(u)f(x-u)du. Nếu trước khi nhân ta không đảo một thừa số thì có tích phân dưới đây, viết là f\star f (dấu sao năm cánh, "pentagram"). Nếu f thực thì f\star f là hàm chẵn, điều không đúng với tích chập nói chung (Bracewell, tr. 40 đến 41).

Chuẩn hóa và cực đại tại gốc

Chia cho giá trị tại gốc ta có \gamma(x) với \gamma(0)=1. Phụ lục chứng minh tự tương quan không vượt giá trị tại gốc: xét \int[f(u)+\varepsilon f(u+x)]^2du\ge0, tam thức bậc hai theo \varepsilon không có nghiệm thực nên b^2-4ac<0, cùng lý luận với bất đẳng thức Schwarz (tr. 41, 48 đến 49). Notebook kiểm trên dãy ngẫu nhiên: tỉ số lớn nhất của các thùy phụ với đỉnh là 0.2821.

Ví dụ 1: f=1-x trên (0,1)

Với 0<x<1: \int f(u)f(u+x)du=\tfrac13-\tfrac x2+\tfrac{x^3}6, giá trị tại gốc là \tfrac13, nên \gamma(x)=1-\tfrac32|x|+\tfrac12|x|^3 và \gamma=0 khi |x|>1. Tại x=0.25: 0.6328; tại x=0.5: 0.3125. Diện tích của tử số là 0.25, đúng bình phương diện tích của f (1/2)^2 (tr. 41 đến 42).

Ví dụ 2: f=e^{-x}H(x)

Với f(x)=e^{-x}H(x): \int f(u)f(u+x)du=\tfrac12e^{-|x|} nên \gamma(x)=e^{-|x|} (tr. 42). Tại x=1: tích phân số cho 0.3679, đúng e^{-1}. Tự tương quan là hàm chẵn dù f chỉ nằm ở nửa dương.

Dữ liệu chạy mãi: đoạn hữu hạn và giới hạn C(x)

Hàm chạy vô hạn thường làm tích phân vô hạn không tồn tại. Cách xử lý là xét một đoạn độ dài X và thay giá trị ngoài đoạn bằng 0, đúng như khi tính trên dữ liệu quan sát hữu hạn. Khi X tăng, \gamma có thể tiến về giới hạn C(x) (tr. 42 đến 43). Khi biến là thời gian, ta viết bằng trung bình thời gian \langle V(t)V(t+\tau)\rangle (tr. 44).

Ba sóng hình sin: pha biến mất

Với V=A\sin(\alpha t+\phi)+B\sin(\beta t+\chi)+C\sin(\gamma t+\psi), tự tương quan giới hạn là chồng ba cosin cùng tần số nhưng biên độ tương đối khác, không còn dấu vết của pha (tr. 43 đến 45). Số đo: A,B,C=2,1,0.5 ở 3, 5, 8 Hz; C(0.05) = 0.4093 với cả hai bộ pha khác nhau (lệch tối đa 2.2e-16).

Tự tương quan không khả nghịch

Vì pha bị mất, không thể từ tự tương quan quay về hàm gốc: tự tương quan là một mất mát thông tin. Trong giao thoa vô tuyến và nhiễu xạ tia X, quan sát tự tương quan dễ hơn quan sát chính hàm, và người ta dồn nhiều công sức để khôi phục thông tin mất (tr. 45). Từ mọi hàm có cùng tự tương quan, hàm gồm toàn cosin có tính duy nhất riêng.

Tương quan bậc ba: giữ chiều của thời gian

Tự tương quan không biết chiều thời gian: nhiễu qua bộ lọc RC e^{-t/RC}H(t) thống kê có thể khác khi phát ngược, nhưng tự tương quan không thấy. Tương quan bậc ba U(x_1,x_2)=\int f(x)f(x+x_1)f(x+x_2)dx thì thấy (tr. 45 đến 46). Bracewell ghi rằng \{1,2,3\} và \{3,2,1\} cho hai mảng khác nhau. Notebook xác nhận: hai dãy có cùng tự tương quan (yes) nhưng mảng bậc ba khác nhau (yes); tâm U(0,0)=\sum f^3 = 36 và tổng mọi phần tử (\sum f)^3 = 216.

Mảng của \{1,2,3\} (hàng x_1=-2,\dots,2, cột x_2=-2,\dots,2):

3 6 9 0 0
6 14 22 6 0
9 22 36 14 3
0 6 14 22 6
0 0 3 6 9

Tương quan chéo và ký hiệu sao năm cánh

Tương quan chéo của g và h là g\star h=\int g(u-x)h(u)du: giống tích chập nhưng g chỉ dịch, không đảo. Đọc là "g quét h". Nó không giao hoán, nhưng g\star h(x)=h\star g(-x) (tr. 46). Với hàm phức phải dùng g^*, và đặt sai vị trí dấu sao sẽ cho liên hợp của kết quả chuẩn (tr. 41, 47).

Notebook kiểm hệ thức với dãy rời rạc: đúng (yes), và g\star h\ne h\star g nói chung (yes).

Phổ năng lượng

Bình phương môđun của biến đổi, |F(s)|^2, là phổ năng lượng. Không có tương ứng một-một với f: cần cả pha. Thông tin mất giống hệt thông tin mất khi dùng tự tương quan, và định lý tự tương quan (chương 6) diễn đạt sự tương đương đó (tr. 47). Notebook kiểm bản rời rạc: biến đổi ngược của |X|^2 trùng tự tương quan vòng (lệch 7.1e-15), và hai tín hiệu khác pha có cùng phổ năng lượng (lệch 2.3e-10).

Phổ năng lượng tích lũy và vạch phổ

Vì |F|^2 là mật độ năng lượng theo s, một vạch phổ mang năng lượng hữu hạn làm mật độ vô hạn. Ta dùng phổ tích lũy \int_0^s|F|^2ds: mỗi vạch là một bước nhảy hữu hạn (hình 3.12, tr. 47 đến 48). Số đo: tín hiệu gồm sóng 50 Hz biên độ 3 cộng nhiễu, phần năng lượng tích lũy trước 49 Hz là 0.005 và sau 51 Hz là 0.951: bước nhảy tại vạch 50 Hz là 0.946 tổng năng lượng.

Dãy xung có tự tương quan tốt: \{1100101\}

Bài tập 1(j), (k) và tài liệu Pettit (1967) trong thư mục nhắc tới dãy xung có tính tự tương quan tốt. Với \{1,1,0,0,1,0,1\}, tự tương quan thẳng có đỉnh 4 tại gốc và thùy phụ lớn nhất 1; tự tương quan vòng bằng 2 tại mọi độ dịch khác 0, gần như phẳng. Tính chất này là lý do dãy như vậy dùng cho radar và đồng bộ.

Điều cần mang theo

• Tích chập có bốn cách nhìn (diện tích, làm mượt, chồng chất, gập) và cư xử như phép nhân: giao hoán, kết hợp, phân phối.
• Với diện tích 1, tích chập nhân diện tích và cộng phương sai; lặp lại nhiều lần tiến về Gauss.
• Tích chập nối tiếp là nhân đa thức; kiểm bằng tổng và độ dài; đảo bằng chia dãy hoặc bảng nghịch đảo.
• Tính bằng máy: mã trực tiếp cần l_fl_g phép nhân; FFT phải đệm số 0 để tránh cuộn vòng.
• Tự tương quan có cực đại tại gốc, chẵn với hàm thực, và mất pha; tương quan bậc ba giữ chiều thời gian; phổ năng lượng |F|^2 tương đương tự tương quan.