加性零均值高斯噪声的产生与信噪比SNR的计算

噪声：引起图像受一些随机误差的影响的退化，噪声可能与图像内容有关，也可能无关。

白噪声：具有常量的功率谱，在所有频率上出现且强度相同的理想的噪声。

高斯噪声：服从高斯正态分布概率的一种白噪声。

在一维下，高斯正太分布概率的密度函数为：

加性噪声：噪声出现与图像信息本身无关的独立于信号的退化噪声。加噪后图像矩阵J可直接表示为：J=I+v，其中v为噪声矩阵。

加性零均值高斯噪声的产生：

设图像I的灰度值范围为[0,255]，取任意数sigma>0，代表噪声的大小程度。

循环：取每对左右相邻的像素点(i,j)与(i,j+1)【如果(i,j)为最右侧点可以跳过或者只对这个点进行加噪】，分别执行如下操作：

1.产生位于[0,1]的随机数r,phi。

2.计算:

3 .加噪后该像素点对应位置计算如下：

J(i,j)=I(i,j)+z1; J(i,j+1)=I(i,j+1)+z2;

4.取下一对左右相邻的像素点直至遍历结束。

遍历完后，J内的像素值可能存在大于255或者小于0的情况，可遍历每个点进行判断，如果大于255则置为255，小于0则置为0，或对J内的像素点进行归一化，另最大的点按比例调整为255，最小的调整为0。

%% Gaussian noise
J=zeros(m,n,d);
sigma=5*rand(1);
for i=1:m
    for j=1:2:n
        r=rand(1);
        phi=rand(1);
        z1=sigma*cos(2*pi*phi)*sqrt(-2*log(r));
        z2=sigma*sin(2*pi*phi)*sqrt(-2*log(r));
        for k=1:d
            J(i,j,k)=I_od(i,j,k)+z1;%I_id为原始图 
            if j+1<=n
                J(i,j+1,k)=I_od(i,j+1,k)+z2;
            end
        end
    end
end
for i=1:m
    for j=1:n
        for k=1:d
            if J(i,j,k)>255%如果大于255则置为255
                J(i,j,k)=255;
            elseif J(i,j,k)<0
                J(i,j,k)=0;
            end
        end
    end
end
J=uint8(J);
figure,
subplot(121),imshow(I_od);
subplot(122),imshow(J);

其中[m,n,d]为原始图像（未加噪图像）的高、宽和色深，彩色则d=3,否则为1。matlab中log()函数表示数学里的ln，而log10()函数表示以10为底，常数e可通过exp(1)获得，pi在matlab中默认定义为3.1416。rand(1)为返回一个1*1的0-1取值的矩阵，也就是返回一个0-1的随机数。在之前为了区分原图像和对应的灰度图像，用I表示灰度图像，I_od表示原图像：

I=imread('lena.jpg');
I_od=I;
[m,n,d]=size(I_od);
if ndims(I)==3%如果为RGB图像则转为灰度
    I=rgb2gray(I);
end

如果本身就是灰度图像，I和I_od为同样大小。matlab中isrgb已不可用，可以替换成ndims(I)==3或者size(I,3)==3。

也可直接通过J=imnoise(I, 'gaussian', 均值，方差)实现。

运行结果如下：

左侧为原始图像，右侧为加噪图像

根据噪声v，可定义信噪比SNR(Signal-to-Noise Radio)，表示图像（或信号）中噪声占总信息数的比例。

首先计算噪声v的所有平方和：

观察到的噪声后信号所有平方和：

信噪比则为E与F的比值：

信噪比又常用他的对数表示，单位为分贝（dB）:

上面加噪后图像J的信噪比计算代码如下：

V=I_od-J;
E=sum(sum(sum(V.^2)));
F=sum(sum(sum(J.^2)));
SNR=F/E;
SNRdB=10*log10(SNR);

直接用加噪后图像和原图像的差值来获得噪声矩阵V，矩阵的平方用V.^2获得，平方和直接sum(sum(V.^2))，因为可能是彩色图像，V为三阶段矩阵，所以要在多进行一次sum（总共三次sum，可能存在一个函数计算矩阵所有元素之和）。

信噪比越大，代表图像质量越好。

除了加性噪声外，很多情况下噪声的大小和图像本身信息有关。

例如：1.乘性噪声:J=I*v，多出现于电视的光栅退化。

2.量化噪声，多出现在量化级别(即灰度级)不足的情况，伪轮廓。

3.冲击噪声：图像像素被个别噪声破坏，被破坏的像素值与邻域内像素值差别较大，椒盐噪声为饱和的冲击噪声，他使得被破坏的像素的像素值要么为255(白盐)，要么为0（黑椒）。

全部可运行代码：（包括直方图那一篇的代码）

%% hist
I=imread('rgb.jpg');
I_od=I;
[m,n,d]=size(I_od);
if ndims(I)==3%如果为RGB图像则转为灰度
    I=rgb2gray(I);
end
%[h,p]=imhist(I);
h=zeros(256,1);
for i=1:m
    for j=1:n
        gray=I(i,j);
        h(gray+1,1)=h(gray+1,1)+1;
    end
end
h=h/sum(h);

K=5;
hf=zeros(256,1);
for i=1:size(h,1)
    for j=-K:K
        if i+j>=1 && i+j<=255
            hf(i,1)=hf(i,1)+h(i+j,1)/(2*K+1);
        end
    end
end
bar(h);
axis([0,255,0,max(h)*1.1]);
figure;
bar(hf);
axis([0,255,0,max(hf)*1.1]);
%% Entropy
Entropy=0;
for i=1:256
    if h(i)~=0
        Entropy=Entropy-h(i)*log(h(i))/log(2);
    end
end
%% Gaussian noise
J=zeros(m,n,d);
sigma=10*rand(1);
for i=1:m
    for j=1:2:n
        r=rand(1);
        phi=rand(1);
        z1=sigma*cos(2*pi*phi)*sqrt(-2*log(r));
        z2=sigma*sin(2*pi*phi)*sqrt(-2*log(r));
        for k=1:d
            J(i,j,k)=I_od(i,j,k)+z1;%I_id为原始彩色图 
            if j+1<=n
                J(i,j+1,k)=I_od(i,j+1,k)+z2;
            end
        end
    end
end
for i=1:m
    for j=1:n
        for k=1:d
            if J(i,j,k)>255%如果大于255则置为255
                J(i,j,k)=255;
            elseif J(i,j,k)<0
                J(i,j,k)=0;
            end
        end
    end
end
J=uint8(J);
figure,
subplot(121),imshow(I_od);
subplot(122),imshow(J);      


V=I_od-J;
E=sum(sum(sum(V.^2)));
F=sum(sum(sum(J.^2)));
SNR=E/F;
SNRdB=10*log10(SNR);