请教：这段程序有何问题，最后如何得出信噪比值？

linisong

请教：这段程序有何问题，最后如何得出信噪比值？
x=0:1:1600;
n=1:length(x);
a=50*exp(-(n-420).^2/144);
d=70*exp(-(n-1200).^2/324);
g=a+d;
h=40*exp(-(n-800).^2/2);
i1=g+h;
r=randn(1,1601);
r=r-mean(r);
r=sqrt(2.56)*r./sqrt(var(r));
ii=i1+r;
y=medfilt1(ii,4);
subplot(3,1,1);plot(x,ii);
subplot(3,1,2);plot(x,y);axis([0 1600 -50 100]);
xd=wden(ii,'minimaxi','s','one',5,'db5');
subplot(3,1,3);plot(x,xd);axis([0 1600 -50 100]);
[row,col,nchannel]=size(i1);
snr=0;
if nchannel==1%gray image
Ps=sum(sum((i1-mean(mean(i1))).^2));%signal power
Pn=sum(sum((i1-y).^2));%noise power
snr=10*log10(Ps/Pn);
elseif nchannel==3%color image
for i=1:3
Ps=sum(sum((i1(:,:,i)-mean(mean(i1(:,:,i)))).^2));%signal power
Pn=sum(sum((i1(:,:,i)-y(:,:,i)).^2));%noise power
snr=snr+10*log10(Ps/Pn);
end
snr=snr/3;
end

=========eight=====================================
请参阅以下帖子：[原创]关于信噪比与向已知信号添加噪声的整理及个人心得
==================================================

[ 本帖最后由 eight 于 2007-1-11 17:45 编辑 ]