%同态滤波
clear all;clc;
a=imread('image.png');
a=rgb2gray(a);
b=double(a);
[m,n]=size(b);
rL=0.5; %默认低频增益
rH=2.0; %默认高频增益
c=1.5; %控制滤波器函数斜面
d0=1800; %截止频率
b1=log(b+1);%取对数
FI=fft2(b1);%二维傅里叶变换
x=floor(m/2);%floor向下取整
y=floor(n/2);
for u=1:m
for v=1:n
D1(u,v)=((u-x).^2+(v-y).^2);%求点到频率矩形中心的距离
H1(u,v)=(rH-rL).*(1-(exp(c*(-D1(u,v)./(d0^2)))))+rL; % 高斯同态滤波
% H1(u,v)=(rH-rL).*(1/(1+(d0/c*(D2(i,j).*1)^2)))+rL; % 巴特沃斯滤波器
end
end
b2=ifft2(H1.*FI);%傅里叶逆变换
b3=real(exp(b2));%指数变换
figure,
subplot(121);imshow(a);title('原图');
subplot(122);
imshow(b3,[]);title('同态滤波后图像');%显示处理后图像
%绘制滤波函数曲线
d=1:1:30000;%设定滤波函数的横坐标范围
H11=(rH-rL)*(1-exp(-c*(d.^2/d0^2)))+rL;
%H11 =(rH-rL)*(1+(d0./d*c).^2).^-1+rL;
figure,
plot(H11);title('滤波函数H(u,v)');
grid on
