求PSNR

function f  =  WPSNR(A,B,varargin)

 %  This function computes WPSNR (weighted peak signal - to - noise ratio) between
 %  two images. The answer  is   in  decibels (dB).
 % 
 %  Using contrast sensitivity function (CSF) to weight spatial frequency
 %  of error image.
 % 
 %  Using:     WPSNR(A,B)
 % 
 %  Written by Ruizhen Liu, http: // www.assuredigit.com 
 
     if  A  ==  B
       error( ' Images are identical: PSNR has infinite value  lilizong ' )
    end

    max2_A  =  max(max(A));
    max2_B  =  max(max(B));
    min2_A  =  min(min(A));
    min2_B  =  min(min(B));

     if  max2_A  >   1   |  max2_B  >   1   |  min2_A  <   0   |  min2_B  <   0 
       error( ' input matrices must have values in the interval [0,1] ' )
    end

    e  =  A  -  B;
     if  nargin < 3 
        fc  =  csf;     %  filter coefficients of CSF
     else 
        fc  =  varargin {1} ;
    end
    ew  =  filter2(fc, e);         %  filtering error with CSF
    
    decibels  =   20 * log10( 1 / (sqrt(mean(mean(ew. ^ 2 )))));
 %     disp(sprintf( ' WPSNR = +%5.2f dB ' ,decibels))
    f = decibels;

 %============= 
function fc  =  csf()
 %============= 
 %  Program to compute CSF
 %  Compute contrast sensitivity function of HVS
 % 
 %  Output:    fc     ---     filter coefficients of CSF
 % 
 %  Reference:
 %     Makoto Miyahara
 %      " Objective Picture Quality Scale (PQS) for Image Coding " 
 %     IEEE Trans. on Comm., Vol  46 , No. 9 ,  1998 .
 % 
 %  Written by Ruizhen Liu, http: // www.assuredigit.com 
 
     %  compute frequency response matrix
    Fmat  =  csfmat;

     %  Plot frequency response
     % mesh(Fmat); pause

     %  compute  2 - D filter coefficient  using  FSAMP2
    fc  =  fsamp2(Fmat);   
     % mesh(fc)

 %=================== 
function Fmat  =  csfmat()
 %=================== 
 %  Compute CSF frequency response matrix
 % 
 %  Called by function csf.m
 % 
 %  Frequency range: the rang of frequency seems to be:
 %      w  =  pi  =  ( 2 * pi * f) / 60 
 %     f  =   60 * w  /  ( 2 * pi),    about  21.2 
 % 
    min_f  =   - 20 ;
    max_f  =   20 ;
    step_f  =   1 ;
    u  =  min_f:step_f:max_f; 
    v  =  min_f:step_f:max_f;
    n  =  length(u);
    Z  =  zeros(n);
     for  i = 1 :n
         for  j = 1 :n
            Z(i,j) = csffun(u(i),v(j));     %  calling function csffun
        end
    end
    Fmat  =  Z;


 %======================== 
function Sa  =  csffun(u,v)
 %======================== 
 %  Contrast Sensitivity Function  in  spatial frequency
 %  This file compute the spatial frequency weighting of errors
 % 
 %  Reference:
 %     Makoto Miyahara
 %      " Objective Picture Quality Scale (PQS) for Image Coding " 
 %     IEEE Trans. on Comm., Vol  46 , No. 9 ,  1998 .
 % 
 %  Input :      u  ---  horizontal spatial frequencies
 %         v  ---  vertical spatial frequencies
 %         
 %  Output:    frequency response
 % 
 %  Written by Ruizhen Liu, http: // www.assuredigit.com 
 
     %  Compute Sa  --  spatial frequency response
     % syms S w sigma f u v
    sigma  =   2 ;
    f  =  sqrt(u. * u + v. * v);
    w  =   2 * pi * f / 60 ;
    Sw  =   1.5 * exp( - sigma ^ 2 * w ^ 2 / 2 ) - exp( - 2 * sigma ^ 2 * w ^ 2 / 2 );

     %  Modification  in  High frequency
    sita  =  atan(v. / (u + eps));
    bita  =   8 ;
    f0  =   11.13 ;
    w0  =   2 * pi * f0 / 60 ;
    Ow  =  (  1   +  exp(bita * (w - w0))  *  (cos( 2 * sita)) ^ 4 )  /  ( 1 + exp(bita * (w - w0)));

     %  Compute final response
    Sa  =  Sw  *  Ow;