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