求psnr值

function psnr=mypsnr(a,b)  
    x=double(a);
    x1=double(b);
    for i=1:256;
        for j=1:256;
            MYa(i,j)=x(i,j)^2;  %after filtering's single
            MYb(i,j)=x(i,j)-x1(i,j);
        end
    end
    d=0;
    e=0;
    for i=1:256;
        for j=1:256;
            d=d+MYa(i,j);
            e=e+MYb(i,j)*MYb(i,j);
        end
    end
    f=log10(d/e);
    g=10*f ;
    psnr=g;

===========================================

function PSNR = mypsnr2(f1, f2)
%计算两幅图像的峰值信噪比

k = 8; %k为图像是表示地个像素点所用的二进制位数，即位深。
fmax = 2.^k - 1;
a = fmax.^2;
e = double(f1) - double(f2);
[m, n] = size(e);
b = sum(e(:).^2);
PSNR = 10*log(m*n*a/b);

=================================================

function a=mypsnr3(A,B)

if A == B
   error('Images are identical: PSNR has infinite value')
end

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

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

error_diff =A - B;
decibels =20*log10(255/(sqrt(mean(mean(error_diff.^2)))));
% disp(sprintf('PSNR =+%5.2f dB',decibels))
a=decibels;

=======================================================

% Function:计算PSNR
function dPSNR = psnrnew(ImageA,ImageB)
if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))
    error(’ImageA <> ImageB’);
    dPSNR = 0;
    return ;
end
M = size(ImageA,1);
N = size(ImageA,2);
d = 0 ;
for i = 1:M
    for j = 1:N
        d = d + (ImageA(i,j) - ImageB(i,j)).^2 ;
    end
end
dPSNR = -10*log10(d/(255*255*M*N)) ;
return

=========================================================

function PSNRDEMO(A,B)

% PURPOSE: To find the PSNR (peak signal-to-noise ratio) between two
%          intensity images A and B, each having values in the interval
%          [0,1]. The answer is in decibels (dB).
%
%          There is also a provision, in EXAMPLE 3 below, for images
%          stored in the interval [0,255], i.e. 256 gray levels.
%
% SYNOPSIS: PSNR(A,B)
%
% DESCRIPTION: The following is quoted from "Fractal Image Compression",
%              by Yuval Fisher et al.,(Springer Verlag, 1995),
%              section 2.4, "Pixelized Data".
%
%              "...PSNR is used to measure the difference between two
%              images. It is defined as
%
%                           PSNR = 20 * log10(b/rms)
%
%              where b is the largest possible value of the signal
%              (typically 255 or 1), and rms is the root mean square
%              difference between two images. The PSNR is given in
%              decibel units (dB), which measure the ratio of the peak
%              signal and the difference between two images. An increase
%              of 20 dB corresponds to a ten-fold decrease in the rms
%              difference between two images.
%             
%              There are many versions of signal-to-noise ratios, but
%              the PSNR is very common in image processing, probably
%              because it gives better-sounding numbers than other
%              measures."
%
% EXAMPLE 1: load clown
%            A = ind2gray(X,map); % Convert to an intensity image in [0,1].
%            B = 0.95 * A;        % Make B close to, but different from, A.
%            PSNR(A,B)            % ---> "PSNR = +33.49 dB"
%
% EXAMPLE 2: A = rand(256); % A is a random 256 X 256 matrix in [0,1].
%            B = 0.9 * A;   % Make B close to, but different from, A.
%            PSNR(A,B)      % ---> "PSNR = +24.76 dB (approx)"
%
% EXAMPLE 3: For images with 256 gray levels: this function PSNR was
%            originally written for matrix-values between 0 and 1,
%            because of MATLAB's preference for that interval.
%
%            However, suppose the matrix has values in [0,255]. Taking
%            Example 1 above, we could change the image to 256 gray levels.
%        
%            load clown
%            A = ind2gray(X,map); % Convert to intensity image in [0,1]
%            AA = uint8(255*A);   % Change to integers in [0,255]
%            BB = 0.95*AA;        % Make BB close to AA.
%
%            Now we must alter the code for this new case. Comment out the
%            existing program (using %) and uncomment the alternative
%            underneath it.
%
%            PSNR(AA,BB)          % ---> "PSNR = +33.56 dB"
%
%            Note the slightly different result from Example 1, because
%            decimal values were rounded into integers.

if A == B
   error('Images are identical: PSNR has infinite value')
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

error_diff = A - B;
decibels = 20*log10(1/(sqrt(mean(mean(error_diff.^2)))));
disp(sprintf('PSNR = +%5.2f dB',decibels))

% if A == B
%    error('Images are identical: PSNR has infinite value')
% end

% max2_A = max(max(A));
% max2_B = max(max(B));
% min2_A = min(min(A));
% min2_B = min(min(B));
%
% if max2_A > 255 || max2_B > 255 || min2_A < 0 || min2_B < 0
%   error('input matrices must have values in the interval [0,255]')
% end

% error_diff = A - B;
% decibels = 20*log10(255/(sqrt(mean(mean(error_diff.^2)))));
% disp(sprintf('PSNR = +%5.2f dB',decibels))

=========================================================

%  copies and that both that copyright notice and this permission notice
% appear in supporting documentation. This software is provided "as is"
% without express or implied warranty. The authors shall not be held
% liable in any event for incidental or consequential damages in
% connection with, or arising out of, the furnishing, performance, or
% use of this program.
%
% If you use the Checkmark software package for your research, please cite:
%
% Shelby Pereira, Sviatoslav Voloshynovskiy, Maribel Madue駉, St閜hane Marchand-Maillet
% and Thierry Pun, Second generation benchmarking and application oriented evaluation,
% In Information Hiding Workshop, Pittsburgh, PA, USA, April 2001.
%
%  http://cui.unige.ch/~vision/Publications/watermarking_publications.html
%
% 
%
% See the also the "Copyright" file provided in this package for
% copyright information about code used in the Checkmark package.
%
function [PSNR,wPSNR]=psnrMetric(a,b,type)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a:
% b:
% type:
% type=1 -> psnr
% type=2 -> nsG
% type=3 -> sgG
% type=4 -> psnr,nsG
% type=5 -> psnr,sgG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wPSNR=[];
if (a==b)
    out=inf;
else
    if size(a,3)==3
        a=[a(:,:,1) a(:,:,2) a(:,:,3)];
        b=[b(:,:,1) b(:,:,2) b(:,:,3)];
    end   
    if type>=1 & type<=5
        if type==2 | type==4
            statistics='nsG';
        else
            statistics='sgG';
        end
        NVF=nvf(a,statistics,150);
        NVF=NVF/max(NVF(:));
        c=NVF.*(a-b).^2;
    end
    c=(a-b).^2;
    PSNR=10*log10(255^2*prod(size(a))/sum(c(:)));
    c=NVF.*c;
    wPSNR=10*log10(255^2*prod(size(a))/sum(c(:)));       
end  

===============================================================

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')
    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 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;


%===================
function Fmat = csfmat()
%===================
% Compute CSF frequency response matrix
% Calling 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;