32f103数字低通滤波函数_数字图像处理：Two-Dimensional Fourier Transform and Filte

╰+哭是因爲堅強的太久メ 2023-01-07 10:08 118阅读 0赞

摘 要：傅立叶变换在图象处理的应用领域中起着极其重要的作用，通过傅里叶变换进行图象分析、图象增强及图象压缩等来获得需要的信息。本实验采用二维快速傅里叶变换对图像进行频域滤波,共编写了一个2-D FFT程序和一个高斯低通滤波器的程序，对图像进行滤波,并通过原图减去高斯低通滤波后的图像实现锐化并且显示图片在二维和三维当中的幅度谱和相位谱。

**欢迎关注、点赞、收藏、私信、交流】共同学习进步**

**一、实验目的**

该实验的目的是开发一个 2-D FFT 程序包。要求程序能完成下面的功能:

(1.a) 用因子乘以输入图像，以实现滤波的中心化变换；

(1.b) 计算傅里叶变换；

(1.c) 用一个实矩阵乘以一个复数矩阵，即用实矩阵中的元素同时乘以复数矩阵对应位置上的复数的实部与虚部。可以通过调用两个图像的乘法程序来实现对应元素的相乘；

(1.d) 计算反傅里叶变换；

(1.e) 结果乘以，并取其实部；

(1.f) 计算频谱。

(1.g) 计算图 4.41(a)的中心化后的傅里叶谱；

(1.f) 使用二维图像和三维图形两种形式显示计算的幅度谱和相位谱；该实验主要实现 4.7.3 节的算法。若试验中使用 MATLAB 语言实现, 则在(1.b)和(1.d)中二维傅里叶变换和反变换程序可以使用 Matlab 的函数 FFT2 和 IFFT2。

(2.a) 编写一个实现公式(4.8-6)所示的高斯低通滤波器的程序。

(2.b) 对图 4.41(a)进行高斯低通滤波以获得图 4.48(c)。

(3.a) 如公式(4.9-1)所示，原始图像减去实验(2.b)中的图像得到锐化后的图像。观察结果图像与图 4.56 中得到的高斯高通滤波结果图的差异， 并解释其中的原因。

(3.b) 调整高斯低通滤波器变量，使(3.a)中的减法锐化图与图 4.56(c)接近为止，并解释你的结果。

**二、技术论述**

1 离散傅里叶变换—DFT

(a)、离散傅里叶变换不是一个新的傅里叶变换形式，它来自于DFS，不过是在时域和频域各取一个周期，再经过周期延拓便可得到整个变换。离散傅里叶变换在时域和频域上都呈现离散的形式，且都是有限长。将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样，而实际上这两组序列都应当被认为是离散周期信号的主值序列。DFT的公式为：

(b)、由于我们要处理的图像是二维的，所以需要用到二维离散傅里叶变换，设f(x,y)是大小为M\*N的数字图像，与一维离散傅里叶变换一样必须在离散变量u和v在u=0,1,2,…,M-1和v=0,1,2,…,N-1范围内求值，二维离散傅里叶变换的公式为：

二维傅里叶反变换公式：

通常二维傅里叶变换的是复函数，因此可以用极坐标的形式来表示，令R和I是F(u,v)的实部和虚部，则：

幅度谱：

相位谱：φ(u,v)=arctan(I(u,v)/R(u,v))

功率谱：

(c)、快速傅里叶变换FFT不是与傅里叶变换不同的另一种变换，而是为了减少DFT运算次数的一种改进算法。它使其成为若干小点数的组合对变换式进行一次次分解，从而减少运算量。通常FFT以2为基数，长度N=2l，当要变换的序列长度不是2的整数次方时，仍然取2为基数的FFT，然后在末位进行补零，使其长度延长至2的整数次方。FFT具有使用非常方便、运算效率高、程序简单等优点。

2 傅里叶谱的中心化处理

在数字图像处理中通常为了清楚地观察F(u,v)的谱线情况需要将原点移到M\*N频域的中心。令：，便得到：

因此若要将图像谱的原点从(0,0)移到图像的中心点(N/2,N/2)处，只要将乘上因子，然后进行傅立叶变换即可。

4 频率域滤波步骤

1.  大小为M\*N的输入图像f=(x+y)，选择填充参数P>=2M-1和Q>=2N-1;
2.  对f(x+y)添加必要数量的0，形成大小P\*Q的填充后的图像;
3.  用乘以移动到其变换中心；
4.  计算自步骤3的DFT，得到;
5.  生成一个实的、对称的滤波函数P\*Q，中心在（P/2,Q/2）处，用阵列相乘形成H和F的乘积；
6.  得到处理后的图像G（x,y）;
7.  通过从G（x,y）的左上象提取M\*N区域，得到最终的处理结果个g(x,y)。

5 高斯低通滤波器

低通滤波是要保留图像中的低频分量而去除图像中的高频分量。图像中的边缘和噪声都对应图像频域中的高频部分，所以通过在频域中的低通滤波可以除去或消弱噪声影响并模糊边缘轮廓，与时域中的平滑方法效果类似。要实现低通滤波需要选一个合适的模板h(u,v)以得到消弱图像频率域H(u,v)的高频分量。高斯滤波器的二维函数为：

式中是距频率矩形中心的距离。是截止频率。

四、实验结果讨论

通过素材Fig4.41(a).jpg进行二维快速傅里叶变换、图像傅里叶频谱、相位谱的计算及高斯低通滤波器的设计与验证。如图4-1所示是原图像傅里叶变换的幅度谱和相位谱，图4-2是原图像的傅里叶变换三维幅度谱和相位谱，图4-3是fftshift函数中心变换后的二维傅里叶幅度谱和fft\_zxzh中心变换后的二维傅里叶图，图4-4是自编程序fft\_zxzh中心变换图的二维傅里叶幅度谱和相位谱，图4-5是fft\_zxzh中心变换的三维傅里叶幅度谱和相位谱，图4-6是中心化变化图和原始图。（程序见附录）

![791e5eade857626beba1292e5dce0d05.png][]

图4-1 原图像傅里叶变换的幅度谱和相位谱

![6ad6df9663d4aae30610be98583ae143.png][]

图4-2原图像的傅里叶变换三维幅度谱和相位谱

![0cf45d40ef914ff808f457d5c0c961ad.png][]

图4-3是fftshift和fft\_zxzh函数中心变换后幅度谱和二维傅里叶图

![5c08071959a9fcd1240cbef68184b769.png][]

图4-4是自编程序fft\_zxzh中心变换图的二维傅里叶幅度谱和相位谱

![57904794a950b269928b1bc097365f72.png][]

图4-5 fft\_zxzh中心变换的三维傅里叶幅度谱和相位谱

![4fc2522dbec6c1403baaf47ce021a5ec.png][]

图4-6 中心化变化图和原始图

设计高斯低通滤波器程序，通过不同的截止频率D0=30、60、160可以对图4.41(a)进行高斯低通滤波以获得滤波后的图像如图4-7所示。通过原图像减去高斯低通滤波后的图像可以得到锐化后的图像如图4-8所示。

![8fe3311b45b58104d2723b96755d879b.png][]

图4-7 不同截止频率时GLPF滤波后的输出图像

![4855a76acf27f6b1fe56a695007b698a.png][]

图4-8（b） 不同截止频率时锐化后的图像

**五、实验结论**

经过对2-D FFT和高斯低通滤波后的图像生进行对比分析可知：

通过图4-1原图像傅里叶变换的幅度谱和相位谱和图4-2原图像傅里叶变换三维的幅度谱和相位谱的对比分析可知，两者在二维的幅度谱和相位谱上没有太大的区别，但是在三维上的幅度谱和相位谱有明显是区别，经过滤波处理后的幅度谱填满了整个坐标轴，其相位谱也更加的集中在三维坐标的中心，将整个相位图进行了移动。

通过图4-3、图4-4、图4-5和图4-6的对比可以发现，图片乘以因子后，整个图片灰度变暗，是因为乘以该因子后，图片的灰度级有一半变为负数，而图片的灰度范围只有\[0 255\],因此灰度级为负的话默认为0，图片再进行取实和反傅里叶变化后，整体图片与原图相比变得比较模糊，那是因为该图片只取了实部的原因。对比前后的频谱图还可以发现，经过滤波后图像的频谱更加的集中了，中间的光点变小，说明达到了滤波的效果。

图4-7为经过高斯低通滤波后的图像，其中使用了不同的截止频率D0为30，60，160，可以发现随着D0的变换，得到的输出图像也会发生变换，当D0变大的时候，输出图像会越来越清晰，当D0=160时，输出图像出现了原图像没有的细节，同时边缘也出现了渐变的黑色区域，因为由于图像是在一个填充区域上实现滤波，故输出的图像边缘出现了渐变的黑色区域。

图4-8为用原图像减去不同截止频率的滤波图像后的锐化效果。可以看到在取较小的值的时候锐化效果更为明显，当=30时，图像要比其他的图片要清晰。观察结果图像与课本图4.56中得到的高斯高通滤波结果图的差异，在取值相同的情况下，原图像减去滤波图像后的输出图像的边缘亮度要高于直接使用高斯滤波器的滤波图像。调整高斯低通滤波器变量时得到的减法锐化图像与图4.56(c)效果接近。

**附录：程序清单**

%1 % 2-D FFT 程序包
    function twoDFFT();
    ima = imread('Fig4.41(a).jpg');
    imshow(ima);
    title('原图像');
    ima_f = fft2(ima); %求原图像二维快速傅里叶变换
    mag_ima = abs(ima_f);          % 在求幅度谱的时候需对傅里叶变换求模
    mag_ima_lg = uint8(log2((mag_ima)))      % 在求幅度谱的时候需对FFT后得到的模取对数
    pha_ima = atan2(imag(ima_f),real(ima_f));  
    % 在求相位谱的时候，需要=使用一个四象限反正切atan2函数来计算
    figure,subplot(1,2,1),imshow(mag_ima_lg,[]);  
    % “[]”是为了使mag_ima_lg的最大最小值分别作为纯白255和纯黑0，中间值映射到0到255之间
    title('原图像二维傅里叶幅度谱');
    subplot(1,2,2),imshow(pha_ima,[]);
    title('原图像二维傅里叶相位谱');
    threeD_mag_ima_lg = log2(mag_ima)+1;
    [x1,y1] = size(threeD_mag_ima_lg);
    [a1,b1]=meshgrid(1:x1,1:y1);  % 先进行meshgrid才能使用mesh函数绘制三维傅立叶图
    figure,subplot(1,2,1),mesh(a1,b1,threeD_mag_ima_lg);   
    shading interp;          % 在网格片内采用颜色插值处理可以使三维图表面更光滑
    title('原像的三维傅里叶幅度谱');
    [x2,y2] = size(pha_ima);
    [a2,b2] = meshgrid(1:x2,1:y2); 
    subplot(1,2,2),mesh(a2,b2,pha_ima); 
    shading interp; 
    title('原图像的三维傅里叶相位谱');
    s = fftshift(ima_f);      % 调用fftshift函数实现图像的中心化变换
    mag_s = abs(s);         % 在求幅度谱的时候需对傅里叶变换求模
    mag_s_lg = uint8(log2(mag_s));  % 在求幅度谱的时候需对FFT后得到的模取对数
    pha_s = atan2(imag(s),real(s));  % 在求相位谱的时候，需要使用一个四象限反正切atan2函数来计算
    figure,subplot(1,2,1);imshow(mag_s_lg,[]);  % “[]”是为了使mag_ima_lg的最大最小值分别作为纯白255和纯黑0，中间值映射到0到255之间
    title('fftshift函数中心变换后的二维傅里叶幅度谱');
    imt = fft_zxzh(ima);   % 调用自编函数fft_zxzh对图像进行中心化变换
    imt_f = fft2(imt);     % 对中心化变换后的图像求二维傅里叶变换
    subplot(1,2,2);imshow(imt_f,[]);
    title('fft_zxzh中心变换后的二维傅里叶图谱');
     
    mag_imt = abs(imt_f);             % 在求幅度谱的时候需对傅里叶变换求模
    mag_imt_lg = uint8(log2(mag_imt)); % 在求幅度谱的时候需对FFT后得到的模取对数
    pha_imt = atan2(imag(imt_f),real(imt_f));  % 在求相位谱的时候，需要=使用一个四象限反正切atan2函数来计算
    figure,subplot(1,2,1),imshow(mag_imt_lg,[]);
    title('fft_zxzh中心变换图的二维傅里叶幅度谱');
    subplot(1,2,2),imshow(pha_imt,[]),
    title('fft_zxzh中心变换图的二维傅里叶相位谱');
     
    threeD_mag_imt_lg = log10(mag_imt)+1;
    [x3,y3] = size(threeD_mag_imt_lg);
    [a3,b3]=meshgrid(1:x3,1:y3);  % 使用函数mesh绘制三维傅立叶图需要先进行meshgrid
    figure,subplot(1,2,1),mesh(a3,b3,threeD_mag_imt_lg); 
    shading interp; 
    %figure,plot3((1:x),(1:y),ak);   shading interp; 
    title('fft_zxzh中心变换图的三维傅里叶幅度谱');
     
    [x4,y4] = size(pha_imt);
    [a4,b4] = meshgrid(1:x4,1:y4);   % 先进行meshgrid才能使用mesh函数绘制三维傅立叶图
    subplot(1,2,2),mesh(a4,b4,pha_imt);  
    shading interp; 
    title('fft_zxzh中心变换图的三维傅里叶相位谱');
    imt_if = ifft2(imt_f);       
    imt_if_abs = abs(imt_if); 
    figure,subplot(1,2,1),imshow(imt_if_abs,[]);
    title('中心化变化图像');
    ima = imread('Fig4.41(a).jpg');
    subplot(1,2,2),imshow(ima);
    title('原始图像');
    
    %2 % fft_zxzh图像中心化转换函数
    function imt = fft_zxzh(ima)    
    [x,y] = size(ima);
    for i = 1:x
        for j = 1:y
          imt(i,j) = ima(i,j)*power((-1),(i+j));%输入图像每个像素点乘以(-1)^x+y因子
        end
    end
    return
    
    %3 %编写一个高斯低通滤波函数对图像进行滤波
    function im=gauslow(ima,D0) 
    D0=160;
    [m,n]=size(ima);
    r=2*m;
    c=2*n;
    %填充图像
    imb=zeros(r,c);
    for i=1:m
        for j=1:n
            imb(i,j)=ima(i,j);
        end
    end
    %利用fftshift函数对频谱进行中心化
    %进行高斯滤波
    imF=fftshift(fft2(imb));
    for i=1:r
        for j=1:c
            H(i,j)=1*exp(-((i-m)*(i-m)+(j-n)*(j-n))/(2*D0*D0));
            imG(i,j)=H(i,j)*imF(i,j);
        end
    end
    img=real(ifft2(imG));
    for i=1:r
        for j=1:c
    img(i,j)=power(-1,i+j)*real(img(i,j));
        end
    end
    img=uint8(img);
    im=zeros(m,n);
    for i=1:m
        for j=1:n
            im(i,j)=img(i,j);
        end
    end
    im=uint8(im); 
    return;
    
    %4 %图像锐化程序
    imb=gauslow(ima,30); %调用高斯低通滤波函数
    imc=gauslow(ima,60);
    imd=gauslow(ima,160);
    figure;
    subplot(1,3,1);imshow(imb);title('D0=30');
    subplot(1,3,2);imshow(imc);title('D0=60');
    subplot(1,3,3);imshow(imd);title('D0=160');
     
    %利用原图减高斯低通滤波后的图像，得到锐化图像
    imb=imresize(imb,[512,512]);
    imb=double(imb);
    imc=imresize(imc,[512,512]);
    imc=double(imc);
    imd=imresize(imd,[512,512]);
    imd=double(imd);
    imB=imsubtract(ima,imb);
    imC=imsubtract(ima,imc);
    imD=imsubtract(ima,imd);
    figure;
    subplot(1,3,1);imshow(imB);title('D0=30时锐化后的图像');
    subplot(1,3,2);imshow(imC);title('D0=60时锐化后的图像');
    subplot(1,3,3);imshow(imD);title('D0=160时锐化后的图像');

[791e5eade857626beba1292e5dce0d05.png]: /images/20221119/ecd1962e00ca4bfd8f74507b29405403.png
[6ad6df9663d4aae30610be98583ae143.png]: /images/20221119/b0def0d78b704cdbaa1cc8018385fe7c.png
[0cf45d40ef914ff808f457d5c0c961ad.png]: https://img-blog.csdnimg.cn/img_convert/0cf45d40ef914ff808f457d5c0c961ad.png
[5c08071959a9fcd1240cbef68184b769.png]: https://img-blog.csdnimg.cn/img_convert/5c08071959a9fcd1240cbef68184b769.png
[57904794a950b269928b1bc097365f72.png]: https://img-blog.csdnimg.cn/img_convert/57904794a950b269928b1bc097365f72.png
[4fc2522dbec6c1403baaf47ce021a5ec.png]: https://img-blog.csdnimg.cn/img_convert/4fc2522dbec6c1403baaf47ce021a5ec.png
[8fe3311b45b58104d2723b96755d879b.png]: https://img-blog.csdnimg.cn/img_convert/8fe3311b45b58104d2723b96755d879b.png
[4855a76acf27f6b1fe56a695007b698a.png]: https://img-blog.csdnimg.cn/img_convert/4855a76acf27f6b1fe56a695007b698a.png