求PSNR

functionf=WPSNR(A,B,varargin)

%ThisfunctioncomputesWPSNR(weightedpeaksignal-to-noiseratio)between
%twoimages.Theanswerisindecibels(dB).
%
%Usingcontrastsensitivityfunction(CSF)toweightspatialfrequency
%oferrorimage.
%
%Using:WPSNR(A,B)
%
%WrittenbyRuizhenLiu,http://www.assuredigit.com

ifA==B
error('Imagesareidentical:PSNRhasinfinitevaluelilizong')
end

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

ifmax2_A>1|max2_B>1|min2_A<0|min2_B<0
error('inputmatricesmusthavevaluesintheinterval[0,1]')
end

e=A-B;
ifnargin<3
fc=csf;%filtercoefficientsofCSF
else
fc=varargin...{1};
end
ew=filter2(fc,e);%filteringerrorwithCSF

decibels=20*log10(1/(sqrt(mean(mean(ew.^2)))));
%disp(sprintf('WPSNR=+%5.2fdB',decibels))
f=decibels;

%=============
functionfc=csf()
%=============
%ProgramtocomputeCSF
%ComputecontrastsensitivityfunctionofHVS
%
%Output:fc---filtercoefficientsofCSF
%
%Reference:
%MakotoMiyahara
%"ObjectivePictureQualityScale(PQS)forImageCoding"
%IEEETrans.onComm.,Vol46,No.9,1998.
%
%WrittenbyRuizhenLiu,http://www.assuredigit.com

%computefrequencyresponsematrix
Fmat=csfmat;

%Plotfrequencyresponse
%mesh(Fmat);pause

%compute2-DfiltercoefficientusingFSAMP2
fc=fsamp2(Fmat);
%mesh(fc)

%===================
functionFmat=csfmat()
%===================
%ComputeCSFfrequencyresponsematrix
%
%Calledbyfunctioncsf.m
%
%Frequencyrange:therangoffrequencyseemstobe:
%w=pi=(2*pi*f)/60
%f=60*w/(2*pi),about21.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);
fori=1:n
forj=1:n
Z(i,j)=csffun(u(i),v(j));%callingfunctioncsffun
end
end
Fmat=Z;


%========================
functionSa=csffun(u,v)
%========================
%ContrastSensitivityFunctioninspatialfrequency
%Thisfilecomputethespatialfrequencyweightingoferrors
%
%Reference:
%MakotoMiyahara
%"ObjectivePictureQualityScale(PQS)forImageCoding"
%IEEETrans.onComm.,Vol46,No.9,1998.
%
%Input:u---horizontalspatialfrequencies
%v---verticalspatialfrequencies
%
%Output:frequencyresponse
%
%WrittenbyRuizhenLiu,http://www.assuredigit.com

%ComputeSa--spatialfrequencyresponse
%symsSwsigmafuv
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);

%ModificationinHighfrequency
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)));

%Computefinalresponse
Sa=Sw*Ow;