How to calculate the wavelet frame of Mexican Hat function

146 Views Asked by At

I was confused these days when I was learning wavelet frame through the famous Daubechies's book - Ten lectures on wavelets. Specially, I only obtained parts of her listed frame bounds results in table 3.1, with the proposition 3.3.2 (edition 1). I estimated the results with some numerical numbers. I list my results below (the same results are shadowed) and you can check the attached MATLAB script if you have time. I do not know where I made some mistakes. Thanks a lot.

the proposition the wavelet frame results

close all;
clear;
clc;


a0_array = [2,power(2,1/2),power(2,1/3),power(2,1/4)];
kb0 = 1.75;
kM = 50;
kK = 10;


A_array = zeros(4,7);
B_array = zeros(4,7);
frame_array = zeros(4,7);

for i_a0 = 1:numel(a0_array)
    a0 = a0_array(i_a0);

    % the first part, sum(|phi((a0^m)*w)|^2)
    i_w = 0;
    sum_energy = zeros(1,numel(1:0.001:a0));
    for w = 1:0.001:a0
        sum_m = 0;
        i_w = i_w + 1;
        for m = -kM:kM
            phi_w = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m))^2)*exp(-((w*(a0^m))^2)/2);
            sum_m = sum_m + abs(phi_w)*abs(phi_w);
        end
        sum_energy(i_w) = sum_m;
        disp(['calculating first part: a0=',num2str(a0),', w=',num2str(w),', sum=',num2str(sum_m)]);
    end


    A0_1 = min(sum_energy);
    B0_1 = max(sum_energy);


    % the second part, sum(sqrt(beta1*beta2))
    i_b0 = 0;
    for b0 = 0.25:0.25:kb0
        i_b0 = i_b0 + 1;
        i_k = 0;
        beda_1 = zeros(1,numel(-kK:kK)-1);
        beda_2 = zeros(1,numel(-kK:kK)-1);
        for k = -kK:kK
            if k == 0
                continue;
            end
            i_k = i_k + 1;

            i_w = 0;
            sum_1_pos = zeros(1,numel(1:0.001:a0));
            sum_2_pos = zeros(1,numel(1:0.001:a0));
            for w = 1:0.001:a0
                i_w = i_w + 1;
                i_m = 0;
                mod_1 = zeros(1,numel(-kM:kM));
                mod_2 = zeros(1,numel(-kM:kM));
                for m = -kM:kM
                    i_m = i_m + 1;
                    r1 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m))^2)*exp(-((w*(a0^m))^2)/2);
                    r2 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m)+2*k*pi/b0)^2)*exp(-((w*(a0^m)+2*k*pi/b0)^2)/2);
                    r3 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m)-2*k*pi/b0)^2)*exp(-((w*(a0^m)-2*k*pi/b0)^2)/2);
                    mod_1_pos(i_m) = abs(r1)*abs(r2);
                    mod_2_pos(i_m) = abs(r1)*abs(r3);
                end
                sum_1_pos(i_w) = sum(mod_1_pos);
                sum_2_pos(i_w) = sum(mod_2_pos);
            end

            i_w = 0;
            sum_1_neg = zeros(1,numel(1:0.001:a0));
            sum_2_neg = zeros(1,numel(1:0.001:a0));
            for w = -a0:0.001:-1
                i_w = i_w + 1;
                i_m = 0;
                mod_1 = zeros(1,numel(-kM:kM));
                mod_2 = zeros(1,numel(-kM:kM));
                for m = -kM:kM
                    i_m = i_m + 1;
                    r1 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m))^2)*exp(-((w*(a0^m))^2)/2);
                    r2 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m)+2*k*pi/b0)^2)*exp(-((w*(a0^m)+2*k*pi/b0)^2)/2);
                    r3 = 2/sqrt(3)*power(pi,-1/4)*((w*(a0^m)-2*k*pi/b0)^2)*exp(-((w*(a0^m)-2*k*pi/b0)^2)/2);
                    mod_1_neg(i_m) = abs(r1)*abs(r2);
                    mod_2_neg(i_m) = abs(r1)*abs(r3);
                end
                sum_1_neg(i_w) = sum(mod_1_neg);
                sum_2_neg(i_w) = sum(mod_2_neg);
            end

            beda_1(i_k) = max([max(sum_1_pos),max(sum_1_neg)]);
            beda_2(i_k) = max([max(sum_2_pos),max(sum_2_neg)]);
            disp(['calculating second part: a0=',num2str(a0),', b0=',num2str(b0),', k=',num2str(k)]);
        end
        deta = sum(sqrt(beda_1.*beda_2));

        A = 2*pi/b0*(A0_1-deta);
        B = 2*pi/b0*(B0_1+deta);
        A_array(i_a0,i_b0) = A;
        B_array(i_a0,i_b0) = B;
        frame_array(i_a0,i_b0) = B/A;
    end
end