数字图像处理实验十理想、高斯、巴特沃斯低通和高通滤波器

利用理想低通、高斯低通、巴特沃斯低通滤波器进行滤波

利用理想高通、高斯高通、巴特沃斯高通滤波器进行滤波

clear;
clc;
%----------------------------图像的灰度化处理-------------------------------
JPG=imread('1.jpg');
JPG=rgb2gray(JPG);

%-----------------------------添加2%椒盐噪声--------------------------------
J=imnoise(JPG,'salt & pepper',0.02);
figure(1);subplot(221);imshow(J);title('原始图像');

%-------------------------添加椒盐噪声后的频谱图----------------------------
F=fft2(J);          %傅里叶变换
F1=log(abs(F)+1);   %取模并进行缩放
% subplot(132);imshow(F1,[]);title('傅里叶变换频谱图');
Fs=fftshift(F);      %FFT频谱中心化，将频谱图中零频率成分移动至频谱图中心
S=log(abs(Fs)+1);    %取模并进行缩放
% subplot(133);imshow(S,[]);title('频移后的频谱图');

%-----------------------------理想低通滤波器--------------------------------
s=fftshift(fft2(im2double(J)));
[m,n0]=size(s);
m0=round(m/2);
n0=round(n0/2);
d0=50; %将理想低通滤波器的截止频率D0设置为50
for i=1:m %双重for循环计算频率点(i,j)与频域中心的距离D(i,j)=sqrt((i-round(m/2)^2+(j-round(n/2)^2))
    for j=1:n0 
        L=sqrt((i-m0)^2+(j-n0)^2);
        if L<=d0  %根据理想低通滤波器产生公式,当D(i,j)<=D0,置为1
            h=1;
        else
            h=0;        %根据理想低通滤波器产生公式,当D(i,j)>D0,置为0
        end
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
%real函数取元素的实部
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(222),imshow(s,[]);
title('理想低通滤波图像'); 

%-----------------------------高斯低通滤波器--------------------------------
s=fftshift(fft2(im2double(J)));
[m,n]=size(s);
d0=30; %将高斯低通滤波器的截止频率D0设置为30
m0=round(m/2);
n0=round(n/2);
for i=1:m
    for j=1:n
        L=sqrt((i-m0)^2+(j-n0)^2);    %根据高斯低通滤波器公式H(u,v)=e^-[D^2(u,v)/(2*D0^2)] 
        h=exp(-(L*L)/(2*(d0^2))); %exp表示以e为底的指数函数
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(223),imshow(s,[]); 
title('高斯低通滤波图像');

%---------------------------巴特沃斯低通滤波器------------------------------
s=fftshift(fft2(im2double(J)));
[m,n]=size(s);%求二维傅里叶变换后图像大小
n0=2;            %将巴特沃斯低通滤波器的阶数n设置为2
d0=30;          %将巴特沃斯低通滤波器的截止频率D0设置为30
m1=round(m/2);
n1=round(n/2);
for i=1:m      %双重for循环计算频率点(i,j)与频域中心的距离D(i,j)=sqrt((i-round(m/2)^2+(j-round(n/2)^2))
    for j=1:n 
        L=sqrt((i-m1)^2+(j-n1)^2);
        h=1/(1+(L/d0)^(2*n0));%根据巴特沃斯低通滤波器公式为1/(1+[D(i,j)/D0]^2n)
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
%real函数取元素的实部
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(224),imshow(s,[]);
title('巴特沃斯低通滤波图像');

%---------------------------高通滤波器相关算法------------------------------
PNG=imread('2.png');
PNG=rgb2gray(PNG);
figure(2);subplot(221);imshow(PNG);title('原始图像');

%-----------------------------理想高通滤波器--------------------------------
s=fftshift(fft2(im2double(PNG)));
[m,n]=size(s);
m0=round(m/2);
n0=round(n/2);
d0=50; %将理想高通滤波器的截止频率D0设置为50
for i=1:m %双重for循环计算频率点(i,j)与频域中心的距离D(i,j)=sqrt((i-round(m/2)^2+(j-round(n/2)^2))
    for j=1:n
        L=sqrt((i-m0)^2+(j-n0)^2);
        if L<=d0  %根据理想高通滤波器产生公式,当D(i,j)<=D0,置为0
            h=0;
        else
            h=1;        %根据理想高通滤波器产生公式,当D(i,j)>D0,置为1
        end
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
%real函数取元素的实部
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(222),imshow(s,[]);
title('理想高通滤波图像');

%-----------------------------高斯高通滤波器--------------------------------
s=fftshift(fft2(im2double(PNG)));
[m,n]=size(s);
d0=10; %将高斯高通滤波器的截止频率D0设置为10
m0=round(m/2);
n0=round(n/2);
for i=1:m
    for j=1:n
        L=sqrt((i-m0)^2+(j-n0)^2);    %根据高斯高通滤波器公式H(u,v)=e^-[D^2(u,v)/2*D0^2] 
        h=1-(exp(-(L^2)/(2*(d0^2)))); %exp表示以e为底的指数函数
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(223),imshow(s,[]); 
title('高斯高通滤波图像');

%---------------------------巴特沃斯高通滤波器------------------------------
s=fftshift(fft2(im2double(PNG)));
[m,n]=size(s);%求二维傅里叶变换后图像大小
n0=2;            %将巴特沃斯高通滤波器的阶数n设置为2
d0=30;          %将巴特沃斯高通滤波器的截止频率D0设置为30
m1=round(m/2);
n1=round(n/2);
for i=1:m      %双重for循环计算频率点(i,j)与频域中心的距离D(i,j)=sqrt((i-round(m/2)^2+(j-round(n/2)^2))
    for j=1:n 
        L=sqrt((i-m1)^2+(j-n1)^2);
        if L==0
            h=0; 
        else
            h=1/(1+(d0/L)^(2*n0));%根据巴特沃斯高通滤波器公式为1/(1+[D0/D(i,j)]^2n)
        end
        s(i,j)=h*s(i,j);%频域图像乘以滤波器的系数
    end
end
%real函数取元素的实部
s=real(ifft2(ifftshift(s)));%最后进行二维傅里叶反变换转换为时域图像
subplot(224),imshow(s,[]);
title('巴特沃斯高通滤波图像');

低通滤波器效果图

高通滤波器效果图