数字图像处理(8):实现FFT快速算法(C语言)
2024-06-21 10:34:52
文章目录
- 1. 实验内容
- 1.1 使用平台及语言
- 1.2 代码流程
- 1.3 FFT、IFFT
- 2. 实验结果
- 2.1 输入图片及其频谱
- 2.2 进行低频滤波
- 2.3 去除直流分量
- 2.4 低频滤波
- 2.5 高频滤波
- 2.6 进一步的高频率波
- 2.7 更进一步的高频滤波
- 3. 遇到的问题及收获
- 3.1 问题一
- 3.2 问题二
- 3.3 问题三
- 附代码：
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++){fread(ptr, 1, 1, fp);Pic[i][j] = *ptr; // 把图像输入到2维数组中,变成矩阵型式 td[i * width + j].re = *ptr;td[i * width + j].im = 0.0;ptr++;}}fclose(fp);FFT2(td, height, width, fd);for (i = 0; i < height; i++){for (j = 0; j < width; j++){double temp = sqrt((fd[i * width + j].re / (height*width)) * (fd[i * width + j].re / (height*width)) + (fd[i * width + j].im / (height*width)) * (fd[i * width + j].im / (height*width)));if (temp > max)max = temp;if (temp < min)min = temp;}}//显示频谱for (i = 0; i < height; i++){for (j = 0; j < width; j++){double temp = sqrt((fd[i * width + j].re / (height*width)) * (fd[i * width + j].re / (height*width)) + (fd[i * width + j].im / (height*width)) * (fd[i * width + j].im / (height*width)));temp = (temp - min) * 25500.0 / (max - min);Pic[i][j] = (unsigned char)temp;//printf("%f\t", temp);}}//频谱更改位置//二四象限置换位置for (i = 0; i < height / 2;i++){for (j = 0; j < width / 2; j++){unsigned char t = Pic[i][j];Pic[i][j] = Pic[height/2 + i][width/2 + j];Pic[height / 2 + i][width / 2 + j] = t;}}//一三象限置换位置for (i = 0; i < height / 2; i++){for (j = width/2; j < width; j++){unsigned char t = Pic[i][j];Pic[i][j] = Pic[height / 2 + i][j - width / 2];Pic[height / 2 + i][j - width / 2] = t;}}//保存频谱图fp = fopen("pinpu.raw", "wb");for (i = 0; i < height; i++){for (j = 0; j < width; j++){fwrite(&Pic[i][j], 1, 1, fp);}}fclose(fp);//对频谱进行处理for (i = 0; i < height; i++){for (j = 0; j < width; j++){if (LOW_PASS == 1){ //低通，滤掉高频。if (((i<DEGREE) && (j<DEGREE)) || ((i>(width - DEGREE)) && (j<DEGREE)) || (i<DEGREE) && (j>(height - DEGREE))) || ((i>(width - DEGREE)) && (j>(height - DEGREE)))){}else{fd[i * width + j].re = 0;fd[i * width + j].im = 0;}}else{ //高通，滤掉低频。即四个角置0if (((i<DEGREE) && (j<DEGREE)) || ((i>(width - DEGREE)) && (j<DEGREE)) || i<DEGREE) && (j>(height - DEGREE))) || ((i>(width - DEGREE)) && (j>(height - DEGREE)))){if (((i<1) && (j<1)) || ((i>(width - 1)) && (j<1)) || ((i<1) && (j>(height - 1))) || ((i>(width - 1)) && (j>(height - 1)))){}else{fd[i * width + j].re = 0;fd[i * width + j].im = 0;}}}}}//保存处理频谱图for (i = 0; i < height; i++){for (j = 0; j < width; j++){double temp = sqrt((fd[i * width + j].re / (height*width)) * (fd[i * width + j].re / (height*width)) + fd[i * width + j].im / (height*width)) * (fd[i * width + j].im / (height*width)));temp = (temp - min) * 25500.0 / (max - min);Pic[i][j] = (unsigned char)temp;//printf("%f\t", temp);}}//频谱更改位置//二四象限置换位置for (i = 0; i < height / 2; i++){for (j = 0; j < width / 2; j++){unsigned char t = Pic[i][j];Pic[i][j] = Pic[height / 2 + i][width / 2 + j];Pic[height / 2 + i][width / 2 + j] = t;}}//一三象限置换位置for (i = 0; i < height / 2; i++){for (j = width / 2; j < width; j++){unsigned char t = Pic[i][j];Pic[i][j] = Pic[height / 2 + i][j - width / 2];Pic[height / 2 + i][j - width / 2] = t;}}fp = fopen("p_pinpu.raw", "wb");for (i = 0; i < height; i++){for (j = 0; j < width; j++){fwrite(&Pic[i][j], 1, 1, fp);}}fclose(fp);IFFT2(td,height, width, fd)for (i = 0; i < height; i++){for (j = 0; j < width; j++){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)td[i * width + j].re;//printf("%f\t", temp);}}fp = fopen("processed.raw", "wb");for (i = 0; i < height; i++){for (j = 0; j < width; j++){fwrite(&Pic[i][j], 1, 1, fp);}}fclose(fp);system("pause");return 0;
}
Fft_ifft.c
#include <math.h>
#include <malloc.h>
#include <string.h>
#include <stdio.h>
#include "fft_ifft.h"/*快速傅里叶变换
TD为时域值,FD为频域值,power为2的幂数*/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{int count;int i, j, k, bfsize, p;double angle;COMPLEX *W, *X1, *X2, *X;/*计算傅里叶变换点数*/count = 1 << power;/*分配运算所需存储器*/W = (COMPLEX *)malloc(sizeof(COMPLEX)*count / 2);X1 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);X2 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);/*计算加权系数*/for (i = 0; i<count / 2; i++){angle = -i* pi * 2 / count;W[i].re = cos(angle);W[i].im = sin(angle);}/*将时域点写入存储器*/memcpy(X1, TD, sizeof(COMPLEX)*count);/*蝶形运算*/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 ;}/*重新排序*/for (j = 0; j<count; j++){p = 0;for (i = 0; i<power; i++){if (j&(1 << i)) p += 1 << (power - i - 1);}FD[j] = X1[p];}/*释放存储器*/free(W);free(X1);free(X2);
}
/*快速傅里叶反变换,利用快速傅里叶变换
FD为频域值,TD为时域值,power为2的幂数*/
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);
}/*************************************************************************
*
* 函数名称：
* Fourier()
*
* 参数:
* COMPLEX* TD -输入的时域序列
* long lWidth -图象宽度
* long lHeight -图象高度
* COMPLEX* FD -输出的频域序列
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数进行二维快速付立叶变换。
*
************************************************************************/
void FFT2(COMPLEX * TD, long lWidth, long lHeight, COMPLEX * FD)
{COMPLEX *TempT, *TempF;// 循环变量long i;long j;// 进行傅里叶变换的宽度和高度（2的整数次方）long w = 1;long h = 1;int wp = 0;int hp = 0;// 计算进行付立叶变换的宽度和高度（2的整数次方）while (w < lWidth){w *= 2;wp++;}while (h < lHeight){h *= 2;hp++;}// 分配内存TempT = (COMPLEX *)malloc(sizeof(COMPLEX)*h);TempF = (COMPLEX *)malloc(sizeof(COMPLEX)*h);// 对y方向进行快速付立叶变换//rgb/*for (i = 0; i < w * 3; i++){// 抽取数据for (j = 0; j < h; j++)TempT[j] = TD[j * w * 3 + i];//rgb// 一维快速傅立叶变换FFT(TempT, TempF, hp);// 保存变换结果for (j = 0; j < h; j++)TD[j * w * 3 + i] = TempF[j];}*///灰度for (i = 0; i < w; i++){// 抽取数据for (j = 0; j < h; j++){TempT[j] = TD[j * w + i];}// 一维快速傅立叶变换FFT(TempT, TempF, hp);// 保存变换结果for (j = 0; j < h; j++){TD[j * w + i] = TempF[j];}}// 释放内存free(TempT);free(TempF);// 分配内存TempT = (COMPLEX *)malloc(sizeof(COMPLEX)*w);TempF = (COMPLEX *)malloc(sizeof(COMPLEX)*w);// 对x方向进行快速付立叶变换//rgb/*for (i = 0; i < h; i++){for (k = 0; k < 3; k++){// 抽取数据for (j = 0; j < w; j++)TempT[j] = TD[i * w * 3 + j * 3 + k];// 一维快速傅立叶变换FFT(TempT, TempF, wp);// 保存变换结果for (j = 0; j < w; j++)FD[i * w * 3 + j * 3 + k] = TempF[j];}}*///灰度for (i = 0; i < h; i++){// 抽取数据for (j = 0; j < w; j++){TempT[j] = TD[i * w + j];}// 一维快速傅立叶变换FFT(TempT, TempF, wp);// 保存变换结果for (j = 0; j < w; j++){FD[i * w + j] = TempF[j];}}// 释放内存free(TempT);free(TempF);
}
/*************************************************************************
*
* 函数名称：
* IFourier()
*
* 参数:
* LPBYTE TD -返回的空域数据
* long lWidth -空域图象的宽度
* long lHeight -空域图象的高度
* COMPLEX* FD -输入的频域数据
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数进行二维快速付立叶逆变换。
*
************************************************************************/
void IFFT2(COMPLEX *TD, long lWidth, long lHeight, COMPLEX * FD)
{COMPLEX *TempT, *TempF;// 循环变量long i;long j;// 进行傅里叶变换的宽度和高度（2的整数次方）long w = 1;long h = 1;int wp = 0;int hp = 0;// 计算进行傅里叶变换的宽度和高度（2的整数次方）while (w < lWidth){w *= 2;wp++;}while (h < lHeight){h *= 2;hp++;}// 计算图像每行的字节数//long lLineBytes = WIDTHBYTES(lWidth * 24);// 分配内存TempT = (COMPLEX *)malloc(sizeof(COMPLEX)*w);TempF = (COMPLEX *)malloc(sizeof(COMPLEX)*w);// 对x方向进行快速付立叶变换//rgb/*for (i = 0; i < h; i++){for (k = 0; k < 3; k++){// 抽取数据for (j = 0; j < w; j++)TempF[j] = FD[i * w * 3 + j * 3 + k];// 一维快速傅立叶变换IFFT(TempF, TempT, wp);// 保存变换结果for (j = 0; j < w; j++)FD[i * w * 3 + j * 3 + k] = TempT[j];}}*///灰度for (i = 0; i < h; i++){// 抽取数据for (j = 0; j < w; j++)TempF[j] = FD[i * w + j];// 一维快速傅立叶变换IFFT(TempF, TempT, wp);// 保存变换结果for (j = 0; j < w; j++)FD[i * w + j] = TempT[j];}// 释放内存free(TempT);free(TempF);TempT = (COMPLEX *)malloc(sizeof(COMPLEX)*h);TempF = (COMPLEX *)malloc(sizeof(COMPLEX)*h);// 对y方向进行快速付立叶变换//rgb/*for (i = 0; i < w * 3; i++){// 抽取数据for (j = 0; j < h; j++)TempF[j] = FD[j * w * 3 + i];// 一维快速傅立叶变换IFFT(TempF, TempT, hp);// 保存变换结果for (j = 0; j < h; j++)TD[j * w * 3 + i] = TempT[j];}*///灰度for (i = 0; i < w; i++){// 抽取数据for (j = 0; j < h; j++)TempF[j] = FD[j * w + i];// 一维快速傅立叶变换IFFT(TempF, TempT, hp);// 保存变换结果for (j = 0; j < h; j++)TD[j * w + i] = TempT[j];}// 释放内存free(TempT);free(TempF);/*for (i = 0; i < h; i++){for (j = 0; j < w * 3; j++){if (i < lHeight && j < lLineBytes)*(TD + i * lLineBytes + j) = FD[i * w * 3 +
j].re;}}*/
}fft_ifft.h
#pragma once
#ifndef COM_H_INCLUDED
#define COM_H_INCLUDED#define pi (double)3.14159265359/*复数定义*/
typedef struct
{double re;double im;
}COMPLEX;/*复数加运算*/
static COMPLEX Add(COMPLEX c1, COMPLEX c2)
{COMPLEX c;c.re = c1.re + c2.re;c.im = c1.im + c2.im;return c;
}/*复数减运算*/static COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{COMPLEX c;c.re = c1.re - c2.re;c.im = c1.im - c2.im;return c;
}/*复数乘运算*/static COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{COMPLEX c;c.re = c1.re*c2.re - c1.im*c2.im;c.im = c1.re*c2.im + c2.re*c1.im;return c;
}
void FFT(COMPLEX * TD, COMPLEX * FD, int power);
void IFFT(COMPLEX *FD, COMPLEX *TD, int power);
void FFT2(COMPLEX * TD, long lWidth, long lHeight, COMPLEX * FD);
void IFFT2(COMPLEX *TD, long lWidth, long lHeight, COMPLEX * FD);#endif有问题多交流，可留言可发邮件，我的邮箱是zhaodongyu艾特pku(这里换成点)edu.cn。
本工程源码已更新至github，欢迎star，欢迎PR:)
数字图像处理(8):实现FFT快速算法(C语言)相关推荐
- Win8Metro(C#)数字图像处理--2.31灰度拉伸算法
Win8Metro(C#)数字图像处理--2.31灰度拉伸算法 原文:Win8Metro(C#)数字图像处理--2.31灰度拉伸算法 [函数名称] 灰度拉伸函数GrayStretchProces ...
- FFT快速傅里叶变换C语言实现信号处理 对振动信号进行实现时域到频域的转换
FFT快速傅里叶变换C语言实现信号处理 对振动信号进行实现时域到频域的转换,可实现FFT8192个点或改成其他FFT1024.4096等等,可以直接运行,运行结果与matlab运行的一致,写好了注释, ...
- 【数字图像处理】经典空域滤波算法
目录 1. 时域.空域.频域 2. 低通(LPF)和高通(HPF) 3. 空域滤波算法 3.1 均值滤波(Mean Filter) 3.2 高斯滤波(Gaussian Filter) 3.3 双边滤波 ...
- 数字图像处理（十一）白平衡算法
文章目录 前言 一.白平衡算法原理 二.算法具体步骤 三.C++代码 四.实验结果 参考 前言 当一副彩色图像数字化后,在显示时颜色有时会看起来有些不正常.这是因为颜色通道中不同的敏感度.增光因子 ...
- 【数字图像处理】双三次插值及其卷积算法（Bicubic Interpolation）
双三次插值数学原理及其卷积算法原理(Bicubic Interpolation) 本文是维基百科上双三次插值的中文翻译,如有侵权会立即删除.本人刚学图像处理,第一次翻译英文文章,水平有限,如有错误还请 ...
- python数字图像处理-图像噪声与去噪算法
图像噪声 椒盐噪声 概述: 椒盐噪声(salt & pepper noise)是数字图像的一个常见噪声,所谓椒盐,椒就是黑,盐就是白,椒盐噪声就是在图像上随机出现黑色白色的像素.椒盐噪声是一种 ...
- C语言数字图像处理---2.1 二值图像形态学算法
本章介绍由数学形态学衍生的二值图像形态学算法,主要包括形态学膨胀.腐蚀.开运算和闭运算四种常用算法,并以此为基础讲解形态学轮廓提取算法,结合C语言编程实现,通俗易懂,图文并茂. [定义与算法] 数学形 ...
- Win8Metro(C#)数字图像处理--2.32图像曝光算法
[函数名称] 图像曝光函数ExposureProcess(WriteableBitmap src,int exposureValue) [函数代码] /// <summary> // ...
- 《数字图像处理》题库4：简答题
前言 这是我在学习数字图像处理这门课程时,从网络上以及相关书籍中搜集到的一些题目, 这些题目主要是针对期末考试的. 做题之前你需要注意以下几点: 这篇文章整理了第4种题型,即简答题. 如果你需要答案, ...
- 数字图像处理:边缘检测(Edge detection)
转载自:https://zhuanlan.zhihu.com/p/59640437 觉得写得通俗易懂,要是每个人的博客都这么人性化.... 写在前面: 本文篇幅较长,用了大量图与公式帮助大家深入理解各 ...
最新文章
- Java中伪造referer来获取数据
- KDE应用如何在GNOME环境下运行？
- React开发（101）:样式处理
- 加密算法使用（四）：AES的使用
- js中函数参数值传递和引用传递
- 模块化日常：开源库与私有库重名
- paip.netbeans断点调试debugger console输出乱码解决方案
- 关于图书馆占座问题的调查
- 11计算机专业vb试题答案,11高三计算机专业VB试题（三）
- 腾讯云云直播、云点播
- python里面while true是什么意思_Python里while True是什么意思?
- 数据恢复中的“老大难”：文件碎片 | 专家约稿
- MOS管和三极管的工作原理对比
- Qt实现屏幕虚拟软键盘
- win7怎么开启文件共享
- 2021高考济南一中成绩查询,2021年济南重点高中名单及排名,济南高中高考成绩排名榜...
- MATLAB FFT算法的应用
- 《李开复自传——世界因你不同》——试读章节pdf下载
- ubuntu18.04向日蔡远程软件安装失败
- 统计学原理 统计中的几个基本概念