Image Denoising -threshold calculation using modified bivariate method

This function used to calculate the threshold used for image denoising using modified threshold method

function [T,sigma1,sigma2] = model3_soft(CD,CH,CH1)
%function used to calculate the threshold and variance using modified 
%bivariate model
%CD and CH : obtained after first stage wavelet decomposition
%CH1: obtained after second stage 2D-wavelet decomposition
sigman = (median(median(abs(CD))))/0.6745;
sigmay1 = 0;
sigmay2 = 0;
[m,n] = size(CD);
[m1,n1] = size(CH1);
%To calculate sigma1 value
for i =1:m
    for j = 1:n
        sigmay1 = sigmay1+((CH(i,j))^2);
    end
end
sigmay1 = sqrt(2)*sigmay1/(m*n);
sigma1 = sqrt(max((((sigmay1))-((sigman)^2)),0));
if sigma1~=0
    T = sqrt(3)*(sigman^2);
else
    T = max(max(abs(CH)));
end
%To calculate sigma2 value
for i =1:m1
    for j = 1:n1
        sigmay2 = sigmay2+((CH1(i,j))^2);
    end
end
sigmay2 = sqrt(2)*sigmay2/(m1*n1);
sigma2 = sqrt(max((((sigmay2))-((sigman)^2)),0));