求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;