巴特沃斯带阻滤波器

Halcom

1. function H = freqfilter_btw(M,N,D0,W,n)
   
2. % 巴特沃斯带阻滤波器
   
3. % input:
   
4. %     滤波器的类型，‘btw’
   
5. %     M,N:频域滤波器的尺寸
   
6. %     D0：带阻滤波器的截止频率
   
7. %     n:巴特沃斯滤波器的阶数
   
8. % output:
   
9. %       H：M xN的矩阵，表示频域滤波器矩阵，数据类型为double，
   
10. if nargin == 5
    
11.         n=1;
    
12. end
    
13. u = -M/2:M/2-1;
    
14. v = -N/2:N/2-1;
    
15. [U,V] = meshgrid(u,v);
    
16. D = sqrt(U.^2+V.^2);
    
17. 
    
18.     H = 1./(1+(D*W./(D.^2-D0^2)).^(2*n));
    
19. end

复制代码

1. function Z = fftfilt2(X,H)
   
2. % 频域滤波
   
3. % 函数输入：
   
4. %       X：输入的空域图像矩阵，double类型
   
5. %       H，频域滤波器，一般为图像X的2倍时较好
   
6. % 函数输出：
   
7. %       Z：输出的空域图像局长呢，数据类型为double类型
   
8. % 二维傅里叶变换
   
9. F = fft2(X,size(H,1),size(H,2));
   
10. % 傅里叶反变换
    
11. Z = H.*F;
    
12. Z = ifftshift(Z);
    
13. Z = abs(ifft2(Z));
    
14. Z = Z(1:size(X,1),1:size(X,2));

复制代码

主程序为：

1. D0 = 50; % 阻止的频率点与频域中心的距离
   
2. W = 3;  % 带宽
   
3. n = 2;  % 阶次
   
4. im = imread('coloredChips.png');           % 原图像
   
5. R = imnoise(im(:,:,1),'gaussian',0,0.01);  % R + 白噪声
   
6. G = imnoise(im(:,:,2),'gaussian',0,0.01);  % G + 白噪声
   
7. B = imnoise(im(:,:,3),'gaussian',0,0.01);  % B + 白噪声
   
8. im = cat(3,R,G,B);                         % 原图像 + 白噪声
   
9. H = freqfilter_btw(2*size(R,1),2*size(R,2),D0,W,n);
   
10. R1 = fftfilt2(R,H);
    
11. G1 = fftfilt2(G,H);
    
12. B1 = fftfilt2(B,H);
    
13. im1 = cat(3,R1,G1,B1);
    
14. im1 = uint8(im1);
    
15. figure('color',[1,1,1])
    
16. subplot(121),imshow(im,[]); title('原始图像')
    
17. subplot(122),imshow(im1,[]); title('巴特沃斯带阻滤波图像')

复制代码