常见噪声的分类与Matlab实现

1.研究噪声特性的必要性

        本文的内容主要介绍了常见噪声的分类与其特性。将噪声建模,然后用模型去实现各式各样的噪声。

        实际生活中的各种照片的老化,都可以归结为以下老化模型。


     这个模型很简单,也可以直接用以下公式来表达。


在频域内,用以下公式区表示。


     根据以上式子,可以看出,老旧照片的复原,主要分为两个任务,一个是去噪;另一个是去卷积,或者称为逆滤波,也就是将老化滤波器做反处理。

     本文首先由噪声类型与其建模。随后的博文,会介绍几种基础的去噪方法和基础的逆滤波方法。

    

2.噪声的实现

      2.1    评价用图像与其直方图

        

      2.2  高斯噪声

        高斯噪声,也称为正态噪声,其统计特性服从正态分布。一种较为泛用的噪声模型。 
        Matlab的实现较为简单,Matlab已经有一个randn(M,N)的函数,用其可以产生出均值为0、方差为1、尺寸为M X N像素的高斯噪声图像。
        用以下程序就可以产生任意均值和方差的高斯噪声。

[plain]  view plain  copy
  1. a = 0;  
  2. b = 0.08;  
  3. n_gaussian = a + b .* randn(M,N);  

         

        2.3 瑞利噪声

        瑞利噪声相比高斯噪声而言,其形状向右歪斜,这对于拟合某些歪斜直方图噪声很有用。

        瑞利噪声的实现可以借由平均噪声来实现。如下所示。


这里的表示均值为0,方差为1的均匀分布的噪声。Matlab里,使用函数rand(M,N)就可以产生一个均值为0,方差为1的均匀噪声。

[plain]  view plain  copy
  1. a = -0.2;  
  2. b = 0.03;  
  3. n_rayleigh = a + (-b .* log(1 - rand(M,N))).^0.5;  

        

       2.4 伽马噪声

         伽马噪声的分布,服从了伽马曲线的分布。伽马噪声的实现,需要使用b个服从指数分布的噪声叠加而来。指数分布的噪声,可以使用均匀分布来实现。


使用若干个(这里用b表示)均匀分布叠加,就可以得到伽马噪声。


当然,当b=1的时候,就可以得到指数噪声了。

[plain]  view plain  copy
  1. a = 25;  
  2. b = 3;  
  3. n_Erlang = zeros(M,N);   
  4.   
  5. for j=1:b  
  6.     n_Erlang = n_Erlang + (-1/a)*log(1 - rand(M,N));  
  7. end  



         2.5 均匀噪声

             如同前面所示,均匀噪声可以由函数rand(M,N)直接产生。


[plain]  view plain  copy
  1. a = 0;  
  2. b = 0.3;  
  3. n_Uniform = a + (b-a)*rand(M,N);  

         2.6 椒盐噪声

         椒盐噪声也成为双脉冲噪声。在早期的印刷电影胶片上,由于胶片化学性质的不稳定和播放时候的损伤,会使得胶片表面的感光材料和胶片的基底欠落,在播放时候,产生一些或白或黑的损伤。事实上,这也可以归结为特殊的椒盐噪声。

        椒盐噪声的实现,需要一些逻辑判断。这里我们的思路是,产生均匀噪声,然后将超过阈值的点设置为黑点,或白点。当然,如果需要拟合电影胶片的损伤的话,可以选用别的类型噪声去拟合。

       

[plain]  view plain  copy
  1. a = 0.05;  
  2. b = 0.05;  
  3. x = rand(M,N);  
  4.   
  5. g_sp = zeros(M,N);  
  6. g_sp = f;  
  7.   
  8. g_sp(find(x<=a)) = 0;  
  9. g_sp(find(x > a & x<(a+b))) = 1;  



3.总结

     本文,实现的几类较为基本的噪声。并给出了其实现的方法,代码在下面。下一篇博文,会进行几个常用去噪滤波器的比较。

[plain]  view plain  copy
  1. close all;  
  2. clear all;  
  3. clc;  
  4.   
  5. f = imread('./original_pattern.tif');  
  6. f = mat2gray(f,[0 255]);  
  7. [M,N] = size(f);  
  8.   
  9. figure();  
  10. subplot(1,2,1);  
  11. imshow(f,[0 1]);  
  12. xlabel('a).Original image');  
  13.   
  14. subplot(1,2,2);  
  15. x = linspace(-0.2,1.2,358);  
  16. h = hist(f,x)/(M*N);  
  17. Histogram = zeros(358,1);  
  18. for y = 1:256  
  19.     Histogram = Histogram + h(:,y);  
  20. end  
  21. bar(-0.2:1/255:1.2,Histogram);  
  22. axis([-0.2 1.2 0 0.014]),grid;  
  23. xlabel('b).The Histogram of a');  
  24. ylabel('Number of pixels');  
  25. %% ---------------gaussian-------------------  
  26. a = 0;  
  27. b = 0.08;  
  28. n_gaussian = a + b .* randn(M,N);  
  29.   
  30. g_gaussian = f + n_gaussian;   
  31.   
  32. figure();  
  33. subplot(1,2,1);  
  34. imshow(g_gaussian,[0 1]);  
  35. xlabel('a).Ruselt of Gaussian noise');  
  36.   
  37. subplot(1,2,2);  
  38. x = linspace(-0.2,1.2,358);  
  39. h = hist(g_gaussian,x)/(M*N);  
  40. Histogram = zeros(358,1);  
  41. for y = 1:256  
  42.     Histogram = Histogram + h(:,y);  
  43. end  
  44. bar(-0.2:1/255:1.2,Histogram);  
  45. axis([-0.2 1.2 0 0.014]),grid;  
  46. xlabel('b).The Histogram of a');  
  47. ylabel('Number of pixels');  
  48.   
  49. %% ---------------rayleigh-------------------  
  50. a = -0.2;  
  51. b = 0.03;  
  52. n_rayleigh = a + (-b .* log(1 - rand(M,N))).^0.5;  
  53.   
  54. g_rayleigh = f + n_rayleigh;   
  55.   
  56. figure();  
  57. subplot(1,2,1);  
  58. imshow(g_rayleigh,[0 1]);  
  59. xlabel('a).Ruselt of Rayleigh noise');  
  60.   
  61. subplot(1,2,2);  
  62. x = linspace(-0.2,1.2,358);  
  63. h = hist(g_rayleigh,x)/(M*N);  
  64. Histogram = zeros(358,1);  
  65. for y = 1:256  
  66.     Histogram = Histogram + h(:,y);  
  67. end  
  68. bar(-0.2:1/255:1.2,Histogram);  
  69. axis([-0.2 1.2 0 0.014]),grid;  
  70. xlabel('b).The Histogram of a');  
  71. ylabel('Number of pixels');  
  72. %% ---------------Erlang-------------------  
  73. a = 25;  
  74. b = 3;  
  75. n_Erlang = zeros(M,N);   
  76.   
  77. for j=1:b  
  78.     n_Erlang = n_Erlang + (-1/a)*log(1 - rand(M,N));  
  79. end  
  80.   
  81. g_Erlang = f + n_Erlang;   
  82.   
  83. figure();  
  84. subplot(1,2,1);  
  85. imshow(g_Erlang,[0 1]);  
  86. xlabel('a).Ruselt of Erlang noise');  
  87.   
  88. subplot(1,2,2);  
  89. x = linspace(-0.2,1.2,358);  
  90. h = hist(g_Erlang,x)/(M*N);  
  91. Histogram = zeros(358,1);  
  92. for y = 1:256  
  93.     Histogram = Histogram + h(:,y);  
  94. end  
  95. bar(-0.2:1/255:1.2,Histogram);  
  96. axis([-0.2 1.2 0 0.014]),grid;  
  97. xlabel('b).The Histogram of a');  
  98. ylabel('Number of pixels');  
  99.   
  100. %% ---------------Exponential-------------------  
  101. a = 9;  
  102. n_Ex = (-1/a)*log(1 - rand(M,N));   
  103.   
  104. g_Ex = f + n_Ex;  
  105.   
  106. figure();  
  107. subplot(1,2,1);  
  108. imshow(g_Ex,[0 1]);  
  109. xlabel('a).Ruselt of Exponential noise');  
  110.   
  111. subplot(1,2,2);  
  112. x = linspace(-0.2,1.2,358);  
  113. h = hist(g_Ex,x)/(M*N);  
  114. Histogram = zeros(358,1);  
  115. for y = 1:256  
  116.     Histogram = Histogram + h(:,y);  
  117. end  
  118. bar(-0.2:1/255:1.2,Histogram);  
  119. axis([-0.2 1.2 0 0.014]),grid;  
  120. xlabel('b).The Histogram of a');  
  121. ylabel('Number of pixels');  
  122.   
  123. %% ---------------Uniform-------------------  
  124. a = 0;  
  125. b = 0.3;  
  126. n_Uniform = a + (b-a)*rand(M,N);  
  127.   
  128. g_Uniform = f + n_Uniform;  
  129.   
  130. figure();  
  131. subplot(1,2,1);  
  132. imshow(g_Uniform,[0 1]);  
  133. xlabel('a).Ruselt of Uniform noise');  
  134.   
  135. subplot(1,2,2);  
  136. x = linspace(-0.2,1.2,358);  
  137. h = hist(g_Uniform,x)/(M*N);  
  138. Histogram = zeros(358,1);  
  139. for y = 1:256  
  140.     Histogram = Histogram + h(:,y);  
  141. end  
  142. bar(-0.2:1/255:1.2,Histogram);  
  143. axis([-0.2 1.2 0 0.014]),grid;  
  144. xlabel('b).The Histogram of a');  
  145. ylabel('Number of pixels');  
  146.   
  147. %% ---------------Salt & pepper-------------------  
  148. a = 0.05;  
  149. b = 0.05;  
  150. x = rand(M,N);  
  151.   
  152. g_sp = zeros(M,N);  
  153. g_sp = f;  
  154.   
  155. g_sp(find(x<=a)) = 0;  
  156. g_sp(find(x > a & x<(a+b))) = 1;  
  157.   
  158. figure();  
  159. subplot(1,2,1);  
  160. imshow(g_sp,[0 1]);  
  161. xlabel('a).Ruselt of Salt & pepper noise');  
  162.   
  163. subplot(1,2,2);  
  164. x = linspace(-0.2,1.2,358);  
  165. h = hist(g_sp,x)/(M*N);  
  166. Histogram = zeros(358,1);  
  167. for y = 1:256  
  168.     Histogram = Histogram + h(:,y);  
  169. end  
  170. bar(-0.2:1/255:1.2,Histogram);  
  171. axis([-0.2 1.2 0 0.3]),grid;  
  172. xlabel('b).The Histogram of a');  
  173. ylabel('Number of pixels');  

原文发于博客:http://blog.csdn.net/thnh169/ 

  • 10
    点赞
  • 71
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值