Nhờ các bác giúp tính convolution của hai hàm.

Nhờ các bác giúp tính convolution của hai hàm.

Em có 2 hàm f1, f2 nhận giá trị trong khoảng [ -7, 7 ] có pdf (probability density function) như hình vẽ dưới. f1 là một phân bố Gauss, f2 cũng là 1 phân bố Gauss nhưng của giá trị làm tròn (round(x)).



Giờ em cần tính pdf của g = f1 +f2. Em biết pdf của g sẽ bằng f1*f2, nghĩa là convolution của f1 và f2. Nhưng em không biết tính ra giá trị cụ thể như thế nào?

Em có matlab trong tay và biết matlab có hàm conv, nhưng hàm này chỉ như kiểu tính nhân 2 đa thức.

Các bác giúp em nhé. Thanks in advance.

Bạn đặt ký hiệu lộn xộn ghê. Có phải là cho biết PDF(f1) và PDF(f2), tìm PDF(f1+f2) ?

Hình vẽ bên trái của bạn không phải là Gaussian PDF.
Thử đoạn code này xem sao:
clear;
x1 = -7;
x2 = 7;
delta = (x2 - x1) / 10000;
x = [x1elta2];
% Gaussian PDF
sigma = 1;
mu = 0;
xA = x;
Na = length(xA);
xA1 = xA(1);
xA2 = xA(end);
fA = (1 / sigma / sqrt(2*pi)) * exp(-(xA - mu).^2 / (sigma^2) / 2);
% quantized-Gaussian probability function
xB1 = round(xA1);
xB2 = round(xA2);
xB = [xB1B2];
Nb = length(xB);
pB = zeros(1,Nb);
for k = 1 : Na
index = round(xA(k)) - xB1 + 1;
pB(index) = pB(index) + fA(k) * delta;
end
%combined PDF
xC1 = xA1 + xB1;
xC2 = xA2 + xB2;
xC = [xC1eltaC2];
Nc = length(xC);
fC = zeros(1,Nc);
for kA = 1 : Na
for kB = 1 : Nb
index = round((xA(kA) + xB(kB) - xC1) / delta) + 1;
fC(index) = fC(index) + fA(kA) * pB(kB);
end
end
figure;
subplot(311);plot(xA,fA);axis([xC1 xC2 0 max(fA)]);
subplot(312);stem(xB,pB);axis([xC1 xC2 0 max(pB)]);
subplot(313);plot(xC,fC);axis([xC1 xC2 0 max(fC)]);
Được werty98 sửa chữa / chuyển vào 20:37 ngày 28/07/2007

Rất cảm ơn bác werty . Em học IT nhưng đợt này lại đi nhận cái thesis toàn Toán là Toán, mà Toán thì em học lâu quá rồi nên chữ nghĩa rơi rớt gần hết. Vấn đề mà em đang vật lộn là ở trong bài báo này ạ:
http://www.springerlink.com/content/06qn4mgtxkajemud/fulltext.pdf
Còn vấn đề mà em hỏi bác là equation (12).
Em có vài điều muốn hỏi bác.
1. Về f1. Xin lỗi bác em viết nhầm, cái plot bên trái không phải là Gaussian. Nó là equation (2) trong bài báo (trong đó nói đây là generalized gaussian distribution). Đồ thị của f1 mà em vẽ tương đương với alpha = 2.7 và beta =0.88:


clear;
x1 = -7;
x2 = 7;
delta = (x2 - x1) / 10000;
x = [x1elta2];
% Generalized Gaussian PDF
alpha = 2.7;
beta = 0.88;
sigma = 1.64;
xA = x;
Na = length(xA);
xA1 = xA(1);
xA2 = xA(end);
fA = (beta/(2*alpha*gamma(1/beta))) * exp( - ((abs(xA)/alpha).^beta));
2. Về f2. f2 là equation (6) trong bài báo. Vì em không nói rõ trong câu hỏi nên bác nghĩ f2 có cùng các tham số Gauss với f1 (đọc code của bác em đoán thế). Trong test của em, sigma_n của (6) có giá trị là 1.64.

3. Với f1 và f2 như trên bác werty có thể sửa lại đoạn code của bác giúp em được không ạ. Em rất cần vụ này.
Rất cám ơn các bác

Được SeeTrocKD sửa chữa / chuyển vào 13:06 ngày 29/07/2007

định trả lời bài ô số soduko mà lại reply nhầm vào đây! Xin lỗi chủ topic nhé.
Được chicken_008 sửa chữa / chuyển vào 16:33 ngày 29/07/2007

Cách lập trình cũng tương tự thế thôi mà. Bạn cứ giữ nguyên fA, pB, thêm fD (tính theo hàm generalized Gaussian PDF) rồi tính fC theo fD và pB (đơn giản chỉ cần thay fA bằng fD).
Một vấn đề nữa là nếu chỉ tính trong đoạn [-7,7] mà hàm fC có heavy tail thì sẽ có sai số đó.

Cám ơn bác werty. Em sẽ thử theo cách bác nói.

Bác werty và các bác cho em hỏi tíêp về vấn đề này.
Vì sợ sai số gây ra bởi heavy tail, em quyết định cho f1 chạy trên toàn trục (từ -256 đến 256 trong bài của em), còn f2 thì chỉ cho chạy trên [-7, 7] (chiếm hơn 99.7% của distribution).
Vấn đề của em là: em cần tính cdf ( cumulative distribution function) của f1*f2.
Nếu tính như kiểu của bác werty nói, thì em chỉ biết giá trị của f1*f2 ở những giá trị rời rạc (discrete), mà như code của bác werty là tại 10.000 điểm trên trục x. Mà như thế thì em tính cdf như thế nào ạ?
Cám ơn các bác

Nếu bạn làm f1 và f2 có resolution khác nhau thì phải cẩn thận coi chừng nhầm lẫn khi tính toán. Tìm hàm cdf thì chỉ việc lấy tích phân lên thôi.
Anyway, không hiểu bạn viết chương trình nhằm mục đích gì ? Hàm f1 và f2 của bạn tương đối đơn giản, có thể tính tích phân ra công thức toán tường minh được.

Hi bác werty,
em cũng đã tính f1*f2 bằng matlab dựa vào định nghĩa (tính tích phân). Tuy nhiên thời gian khá lâu. Mà phép tính convolution này em sẽ phải chạy vài trăm nghìn lần (có thể là vài triệu lần) nên em cần cách tính nhanh nhanh cho hàm này.
Mô tả ngắn gọn my task: các giá trị alpha, beta là các tham số của nhiễu của một bức ảnh (image noise). Nhiễu này bị tạo ra bởi một steganographic method nào đó. Nhiệm vụ của em là xác định các tham số này khi chỉ biết ảnh nhiễu chứ không biết ảnh gốc. Vì vậy, em cần thực hiện 1 grid search với các giá trị của alpha và beta rồi Chi-square gof test để tìm ra cặp giá trị (alpha, beta) tốt nhất: f1*f2 mà mình tính là expected distribution, mình muốn nó khớp với observed distribution của bức ảnh nhiễu.
Thanks

Được SeeTrocKD sửa chữa / chuyển vào 12:06 ngày 03/08/2007
Được SeeTrocKD sửa chữa / chuyển vào 13:08 ngày 03/08/2007