（4）MATLAB 频率域图像增强

   频域滤波可以用来消除周期噪声。

    傅里叶变换

    图像较平滑，低频部分对应的幅值较大，图像灰度变化越剧烈，其频谱高频分量较强

I1 = imread('cell.tif');  %读入原图像

fcoef = fft2(double(I1)); %做fft变换 
spectrum = fftshift(fcoef); %将零点移到中心
temp =log(1+abs(spectrum)); %对幅值做对数变换以压缩动态范围

subplot(1,2,1);
imshow(temp,[]);
title('FFT');
subplot(1,2,2);
imshow(I1);
title('Source')

  幅度谱反映图像整体上各个频率分量的相对强度。相位谱对应图像内容。

% 读取图片
A = imread('../beauty.jpg');
B = imread('../cat.jpg');
% 求傅立叶变换
Af = fft2(double(A));
Bf = fft2(double(B));
% 分别求幅度谱和相位谱
AfA = abs(Af);
AfB = angle(Af);
BfA = abs(Bf);
BfB = angle(Bf);
% 交换相位谱并重建复数矩阵
AfR = AfA .* cos(BfB) + AfA .* sin(BfB) .* i;
BfR = BfA .* cos(AfB) + BfA .* sin(AfB) .* i;
% 傅立叶反变换
AR = abs(ifft2(AfR));
BR = abs(ifft2(BfR));
% 显示图像
subplot(2,2,1);
imshow(A);
title('美女原图像');
subplot(2,2,2);
imshow(B);
title('猫的原图像');
subplot(2,2,3);
imshow(AR, []);
title('美女的幅度谱和猫的相位谱组合');
subplot(2,2,4);
imshow(BR, []);
title('猫的幅度谱和美女的相位谱组合');

理想低通滤波

function out = imidealflpf(I, freq)
% imidealflpf函数			构造理想的频域低通滤波器
% I参数				输入的灰度图像
% freq参数				低通滤波器的截止频率
% 返回值：out – 指定的理想低通滤波器
[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        if (sqrt(((i-M/2)^2+(j-N/2)^2))>freq)
            out(i,j)=0;
        end
    end
end

I = imread('../baby_noise.bmp'); %读入原图像

% 生成滤镜
ff = imidealflpf(I, 20);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,1);
imshow(I);
title('Source');

% 计算FFT并显示
temp = fft2(double(I));
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,1);
imshow(temp, []);
title('Source');

figure (1);
subplot(2,2,2);
imshow(out);
title('Ideal LPF, freq=20');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,2);
imshow(temp, []);
title(' Ideal LPF, freq=20');

% 生成滤镜
ff = imidealflpf(I, 40);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,3);
imshow(out);
title('Ideal LPF, freq=40');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,3);
imshow(temp, []);
title(' Ideal LPF, freq=40');

% 生成滤镜
ff = imidealflpf(I, 60);
% 应用滤镜
out = imfreqfilt(I, ff);

figure (1);
subplot(2,2,4);
imshow(out);
title('Ideal LPF, freq=60');

% 计算FFT并显示
temp = fft2(out);
temp = fftshift(temp);
temp = log(1 + abs(temp));
figure (2);
subplot(2,2,4);
imshow(temp, []);
title(' Ideal LPF, freq=60');

高斯低通滤波器

function out = imgaussflpf(I, sigma)
% imgaussflpf函数     		构造频域高斯低通滤波器
% I参数				输入的灰度图像
% sigma参数			高斯函数的Sigma参数

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);
    end
end

高斯高通滤波器

function out = imgaussfhpf(I, sigma)
% imgaussfhpf函数			构造频域高斯高通滤波器
% I参数				输入的灰度图像
% sigma参数			高斯函数的Sigma参数

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = 1 - exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);
    end
end

高斯带阻滤波器

function out = imgaussfbrf(I, freq, width)
% imidealflpf函数 		构造频域高斯带阻滤波器
% I参数				输入的灰度图像
% freq参数				阻带中心频率
% width参数			阻带宽度

[M,N] = size(I);
out = ones(M,N);
for i=1:M
    for j=1:N
        out(i,j) = 1-exp(-0.5*((((i-M/2)^2+(j-N/2)^2)-freq^2)/(sqrt(i.^2+j.^2)*width))^2);
    end
end