数字图像处理(8):实现FFT快速算法(C语言)

1. 实验内容

1.1 使用平台及语言

使用平台：VS2015

语言：C语言

1.2 代码流程

其中，对频谱图像进行了一三、二四象限的对调，这样便于观察分析。

1.3 FFT、IFFT

对图像进行操作要进行两次，先对行进行FFT，再对列进行FFT，顺序反过来也可以。

FFT最主要的是蝶形运算

    /*蝶形运算*/
    for (k = 0; k<power; k++)
    {
   
   
        for (j = 0; j<1 << k; j++)
        {
   
   
            bfsize = 1 << (power - k);
            for (i = 0; i<bfsize / 2; i++)
            {
   
   
                p = j*bfsize;
                X2[i + p] = Add(X1[i + p], X1[i + p + bfsize / 2]);
                X2[i + p + bfsize / 2] = Mul(Sub(X1[i + p],X1[i + p + bfsize / 2]), W[i*(1 << k)]);
            }
        }
        X = X1;
        X1 = X2;
        X2 = X ;
    }

IFFT计算过程：

数据取共轭,然后fft,结果再取共轭后除以N

void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
   
   
    int i, count;
    COMPLEX *x;
    /*计算傅里叶反变换点数*/
    count = 1 << power;
    /*分配运算所需存储器*/
    x = (COMPLEX *)malloc(sizeof(COMPLEX)*count);
    /*将频域点写入存储器*/
    memcpy(x, FD, sizeof(COMPLEX)*count);
    /*求频域点的共轭*/
    for (i = 0; i<count; i++)
    {
   
   
        x[i].im = -x[i].im;
    }
    /*调用快速傅里叶变换*/
    FFT(x, TD, power);
    /*求时域点的共轭*/
    for (i = 0; i<count; i++)
    {
   
   
        TD[i].re /= count;
        TD[i].im = -TD[i].im / count;
    }
    /*释放存储器*/
    free(x);
}

2. 实验结果

2.1 输入图片及其频谱

2.2 进行低频滤波

将频谱图中心20*20区域置零，进行IFFT变换，得到结果

分析：得到的结果比较恐怖，和预期不符。怀疑是把直流分量置零的原因。

2.3 去除直流分量

对频谱中心2*2置0（即去除直流成分）结果

2.4 低频滤波

进行低频滤波，将频谱图中心20*20区域置零（除了中间的直流成分），进行IFFT变换，得到结果

分析：这是高通滤波，低频成分被去除。图片的效果很明显，变化不大的地方被抑制，留下的都是变化快的地方、边缘部分（高频部分）。

2.5 高频滤波

进行高频滤波，除了频谱中心300*300区域其他置零，进行IFFT变换，得到结果

分析：这是低通滤波，高频成分被去除。图片里变化快的地方（树的纹理、博雅塔受到影响）。

2.6 进一步的高频率波

上个效果不是很明显，重新进行高频滤波，除了频谱中心200*200区域其他置零，进行IFFT变换，得到结果

分析：对比博雅塔部分——此次的结果和上次没有太大差别，只是云那儿受到影响多了一点。怀疑是加的这点效果对博雅塔没什么效果（因为这儿本来就高频比较高吧）

2.7 更进一步的高频滤波

上个效果加强其实也不是很明显，重新进行高频滤波，除了频谱中心100*100区域其他置零，进行IFFT变换，得到结果

分析：没想到这么丑，感觉这次的结果可以和高频滤波结果对比观看——

很明显二者重点的区域完全相反。

3. 遇到的问题及收获

3.1 问题一

最开始有一次遇到了这个问题（堆栈溢出问题），遇过这个坑，所以解决得比较顺利——

这个问题的定位在于最开始定义了数组

    COMPLEX td[height*width];
    COMPLEX fd[height*width];

其中

typedef struct
{
        double re;
        double im;
}COMPLEX;

图像大小是512*512的，这么一算，这两句需要的堆栈就是
2 * 512 * 512 * 2 * 8Bit(size of double) = 8M

解决方法有两种：一是设置堆栈增大容量；二是开辟内存存放数据。我用的第二种。即

td =(COMPLEX*)malloc(height*width*sizeof(COMPLEX));
if (td == NULL)
    return -1;
fd = (COMPLEX*)malloc(height*width*sizeof(COMPLEX));
if (fd == NULL)
    return -1;

3.2 问题二

关于频谱图归一化问题,在显示频谱图时，尝试进行了归一化至0~255

temp = (temp - min) * 255.0 / (max - min);

但是效果很不好——

如图，只有中间有个白点，别的看不出来。不得已将temp的值乘以100后才能看到比较好的效果。

感觉这一问题应该是老师上课讲的直流分量数值太大的原因。致使即使归一化也还是不能很好的展现高频分量。

3.3 问题三

进行IFFT变换后，保存图像时想当然的将复数的幅度的值给了图像，代码如下——

double temp = sqrt((td[i * width + j].re) * (td[i * width + j].re) + (fd[i * width + j].im) * (fd[i * width + j].im));
Pic[i][j] = (unsigned char)temp;

未对频域进行处理，只对图像进行FFT和IFFT，输出结果如下

更改代码，只将实部赋值给图像

Pic[i][j] = (unsigned char)td[i * width + j].re;

效果正常

附代码：

main.c

#include <stdio.h>  
#include <stdlib.h>  
#include <string.h>  
#include "fft_ifft.h"
#include <math.h>

#define height  512  
#define width   512 
#define LOW_PASS 1   //是否为低通
#define DEGREE   150     //滤波程度
typedef unsigned char  BYTE;    // 定义BYTE类型，占1个字节  

int main(void)
{
   
   
    FILE *fp = NULL;
    //  BYTE Pic[height][width];
    BYTE *ptr;
    BYTE **Pic = new BYTE *[height];
    for (int i = 0; i != height; ++i)
    {
   
   
        Pic[i] = new BYTE[width];
    }
    int i, j;
    double max = 0;
    double min = 255;
    COMPLEX *td = NULL;
    COMPLEX *fd = NULL;
    //COMPLEX td[height*width];
    //COMPLEX fd[height*width];
    td = (COMPLEX*)malloc(height*width*sizeof(COMPLEX));
    if (td == NULL)
        return -1;
    fd = (COMPLEX*)malloc(height*width*sizeof(COMPLEX));
    if (fd == NULL)
        return -1;
    fp = fopen("weiminglake512x512.raw", "rb");
    ptr = (BYTE*)malloc(width * height * sizeof(BYTE));//创建内存
    for (i = 0; i < height; i++)
    {
   
   
        for (j = 0; j < width; j++)
        {