【经典算法实现 45】BMP图像的FFT快速傅里叶变换及 IFFT逆变换（C语言）
在前文《【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》中，
我们实现了二维数组的傅里叶变换，在它的基础上，我们来实现图像的FFT快速傅里叶变换及 IFFT逆变换。
主要实现步骤为：
- 从RGB24色 bmp 图片中读取图片信息，将 RGB 转换为 YUV分量
y = (unsigned char)( ( 66 * r + 129 * g + 25 * b + 128) >> 8) + 16 ;
u = (unsigned char)( ( -38 * r - 74 * g + 112 * b + 128) >> 8) + 128 ;
v = (unsigned char)( ( 112 * r - 94 * g - 18 * b + 128) >> 8) + 128 ;
- 将转换后的Y分量，先经过傅里叶变换，变换后的结果，去中心化后，保存复数的模为一张Y分量图片。
- 将Y分量图片，再经过傅里叶逆变换，将结果再保存为一张Y分量图片，对比变换前后的图片。
本文链接：《【经典算法实现 45】C语言实现BMP图像的FFT快速傅里叶变换及 IFFT逆变换》
本文涉及的C代码及图片素材已上传到CSDN：《bmp图片的快速傅里叶变换.zip》
一、图像的FFT快速傅里叶变换及 IFFT逆变换 代码
代码主要包括，傅里叶变换部分 及 BMP 文件读写部分，完整如下：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <float.h>
// 快速傅里叶变换的数据必须是 2^n
int NUM_x, NUM_y;
// 将傅里叶变换后的四角移到中心
void fft_shift(struct _complex *src)
{
int x, y, a, b;
struct _complex tmp;
for(y = 0; y<NUM_y/2; y++){
for(x = 0; x<NUM_x; x++){
a = y * NUM_x + x;
b = ( (y + NUM_y/2)%NUM_y )*NUM_x + (NUM_x/2 + x)%NUM_x;
tmp = src[a];
src[a] = src[b];
src[b] = tmp;
}
}
}
// 以 {0，2，4，6， 1，3，5，7} 排序，
// 数组传参使用 &f[y][x] 的形式传入，取数据时使用 f[y*NUM_x + x].x 形式获取
// split_array(&f[u][0], &F[u][0], NUM_x, 0, 0); // 对 x分组
// split_array(&f[u][0], &F[u][0], 0, NUM_y, 1); // 对 y分组
// flag: 0 时，对x 分组
// flag: 1 时，对y 分组
void split_array(struct _complex *src, struct _complex *dst , int x_n, int y_n, int flag)
{
int i;
struct _complex t[flag == 0 ? x_n/2 : y_n/2], *s = src, *d = dst;
if(flag == 0){
if(x_n <= 1) return;
for(i = 0; i<x_n/2 ; i++){
t[i].x = s[i*2 + 1].x; // 暂存奇数项
t[i].y = s[i*2 + 1].y;
d[i].x = s[i*2].x; // 拷贝偶数项到低位
d[i].y = s[i*2].y;
}
for(i = 0; i<x_n/2 ; i++){
d[i + x_n/2].x = t[i].x; // 拷贝奇数项到高位
d[i + x_n/2].y = t[i].y;
}
split_array(dst, dst, x_n/2, y_n, flag);
split_array(dst+x_n/2, dst+x_n/2, x_n/2, y_n, flag);
}
else
{
if(y_n <= 1) return;
for(i = 0; i<y_n/2 ; i++){
t[i].x = s[(i*2 + 1)*NUM_x ].x; // 暂存奇数项
t[i].y = s[(i*2 + 1)*NUM_x ].y;
d[i*NUM_x ].x = s[(i*2)*NUM_x ].x; // 拷贝偶数项到低位
d[i*NUM_x ].y = s[(i*2)*NUM_x ].y;
}
for(i = 0; i<y_n/2 ; i++){
d[(i+y_n/2)*NUM_x ].x = t[i].x; // 拷贝奇数项到高位
d[(i+y_n/2)*NUM_x ].y = t[i].y;
}
split_array(dst, dst, x_n, y_n/2, flag);
split_array(dst+NUM_x*y_n/2, dst+NUM_x*y_n/2, x_n, y_n/2, flag);
}
}
void printf_array(struct _complex *dst)
{
int x,y;
for(y = 0; y<NUM_y; y++)
{
for(x = 0; x<NUM_x; x++){
printf("[%7.2f+%7.2fj] ", dst[y*NUM_x +x].x, dst[y*NUM_x +x].y);
}
printf("\n");
}
}
// src:源数组 &f[y][x] dst:结果数组&F[y][x] flag: 1:正变换 -1 逆变换
void fft(struct _complex *src, struct _complex *dst, int flag)
{
int y, x, i, u, k , n;
double wu;
struct _complex w, a0, a1, t;
clock_t start, end;
start = clock();
// 傅里叶正变换时，将中心平移到四边
if(flag == -1)
fft_shift(src);
///
// 对x每一行做傅里叶变换
for(y = 0; y<NUM_y; y++){
// 先分割数组
split_array(&src[y*NUM_x+0], &dst[y*NUM_x+0], NUM_x, 0, 0); // 对 x分组
// 对 f[y][] 这一组数进行傅里叶变换
for(i = 0; i<log2(NUM_x); i++){ //计算次数为 2^n = num，即n = log2^num
// 每次计算的间隔是 2^n,分别为 1，2，4，8
n = 2 * pow(2, i); // 本轮一组个数为 2 * 2^n，分别为 2,4,8,则好3轮
for(k = 0; k<NUM_x/n; k++){ // num/n 为当前的组数，分别为 4，2，1
for(u=0; u<n/2; u++){ // 对每组进行计算, a0 和 b0 的个数分别为 n/2
wu = -1 * 2 * M_PI * u / n; // 计算旋转因子
w.x = cos(wu);
w.y = flag * sin(wu); // 如果是傅里叶逆变换，此处 flag = -1
a0 = dst[y*NUM_x + k*n + u]; // 奇数项 [y][k*n+u]
a1 = dst[y*NUM_x + k*n + u + n/2]; // 偶数项 [y][k*n+u+n/2]
t.x = w.x*a1.x - w.y*a1.y;
t.y = w.x*a1.y + w.y*a1.x;
dst[y*NUM_x + k*n + u].x = a0.x + t.x; // F[u] = A0 + wA1
dst[y*NUM_x + k*n + u].y = a0.y + t.y;
dst[y*NUM_x + k*n + u + n/2].x = a0.x - t.x; // F[u+n/2] = A0 - wA1
dst[y*NUM_x + k*n + u + n/2].y = a0.y - t.y;
}
}
}
if(flag == 1) // 正向时，除一个X长度，避免数据过大，本来不应该除的
for(u = 0; u<NUM_x; u++){
dst[y*NUM_x + u].x /= NUM_x;
dst[y*NUM_x + u].y /= NUM_x;
}
}
/*
// 打印当变换结果
end = clock();
printf("\n\n每行傅里叶%s变换结果为： 耗时%fs, NUM=%d x %d\n",flag==1?"":"逆", (double)(end-start)/CLOCKS_PER_SEC, NUM_x, NUM_y);
printf_array(dst);
*/
///
// 对y每一列做傅里叶变换
for(x = 0; x<NUM_x; x++){
// 先分割数组
split_array(&dst[0*NUM_x+x], &dst[0*NUM_x+x], 0, NUM_y, 1); // 对 y分组
// 对 f[][x] 这一组数进行傅里叶变换
for(i = 0; i<log2(NUM_y); i++){ //计算次数为 2^n = num，即n = log2^num
// 每次计算的间隔是 2^n,分别为 1，2，4，8
n = 2 * pow(2, i); // 本轮一组个数为 2 * 2^n，分别为 2,4,8,则好3轮
for(k = 0; k<NUM_y/n; k++){ // num/n 为当前的组数，分别为 4，2，1
for(u=0; u<n/2; u++){ // 对每组进行计算, a0 和 b0 的个数分别为 n/2
wu = -1 * 2 * M_PI * u / n; // 计算旋转因子
w.x = cos(wu);
w.y = flag * sin(wu); // 如果是傅里叶逆变换，此处 flag = -1
a0 = dst[(k*n + u)*NUM_x + x ]; // 奇数项 [k*n+u][x]
a1 = dst[(k*n + u + n/2)*NUM_x + x ]; // 偶数项 [(k*n + u + n/2)*NUM_y + x ][x]
t.x = w.x*a1.x - w.y*a1.y;
t.y = w.x*a1.y + w.y*a1.x;
dst[(k*n + u)*NUM_x + x ].x = a0.x + t.x; // F[u] = A0 + wA1
dst[(k*n + u)*NUM_x + x ].y = a0.y + t.y;
dst[(k*n + u + n/2)*NUM_x + x ].x = a0.x - t.x; // F[u+n/2] = A0 - wA1
dst[(k*n + u + n/2)*NUM_x + x ].y = a0.y - t.y;
}
}
}
if(flag == -1)
for(u = 0; u<NUM_y; u++){
dst[u*NUM_x + x].x /= NUM_y;
dst[u*NUM_x + x].y /= NUM_y;
}
}
// 傅里叶正变换时，将四边平移到中心
if(flag == 1)
fft_shift(dst);
// 打印当变换结果
end = clock();
printf("\n傅里叶%s变换结束，耗时%fs, NUM=%d x %d\n",flag==1?"":"逆", (double)(end-start)/CLOCKS_PER_SEC, NUM_x, NUM_y);
//printf_array(dst);
}
/
// bmp代码部分
/
typedef struct bmp_head
{
//__int16 bfType; //2Bytes，必须为"BM"，即0x424D 才是Windows位图文件
__int32 bfSize; //4Bytes，整个BMP文件的大小
__int16 bfReserved1; //2Bytes，保留，为0
__int16 bfReserved2; //2Bytes，保留，为0
__int32 bfOffBits; //4Bytes，文件起始位置到图像像素数据的字节偏移量
} bmp_head;
typedef struct bmp_info
{
__int32 biSize; //4Bytes，INFOHEADER结构体大小，存在其他版本I NFOHEADER，用作区分
__int32 biWidth; //4Bytes，图像宽度（以像素为单位）
__int32 biHeight; //4Bytes，图像高度，+：图像存储顺序为Bottom2Top，-：Top2Bottom
__int16 biPlanes; //2Bytes，图像数据平面，BMP存储RGB数据，因此总为1
__int16 biBitCount; //2Bytes，图像像素位数
__int32 biCompression; //4Bytes，0：不压缩，1：RLE8，2：RLE4
__int32 biSizeImage; //4Bytes，4字节对齐的图像数据大小
__int32 biXPelsPerMeter; //4 Bytes，用象素/米表示的水平分辨率
__int32 biYPelsPerMeter; //4 Bytes，用象素/米表示的垂直分辨率
__int32 biClrUsed; //4 Bytes，实际使用的调色板索引数，0：使用所有的调色板索引
__int32 biClrImportant; //4 Bytes，重要的调色板索引数，0：所有的调色板索引都重要
}bmp_info;
typedef struct bmp_color
{
char rgbBlue; //指定蓝色强度
char rgbGreen; //指定绿色强度
char rgbRed; //指定红色强度
char rgbReserved; //保留，设置为0
} bmp_color;
unsigned char bfType[2];
bmp_head bmp_h;
bmp_info bmp_i;
unsigned char *bmp_data;
struct _complex *bmp_data_F;
void open_bmp(char * url)
{
int i, j , r, g ,b, y;
FILE *fp = fopen(url, "rb+");
char str[10];
printf("\n\n文件《%s》信息：\n", url);
//读取 bfType
fread(&bfType, 1, 2, fp);
printf("[bfType]: %c%c\n", bfType[0], bfType[1]);
//读取 bmp_head 信息
fread(&bmp_h, sizeof(bmp_head), 1, fp);
printf("[bfSize]: %d\n[bfOffBits]:%d\n",bmp_h.bfSize, bmp_h.bfOffBits);
bmp_h.bfReserved1 = 0;
bmp_h.bfReserved2 = 0;
//读取 bmp_info 信息
fread(&bmp_i, sizeof(bmp_info), 1, fp);
printf("\n[biSize]: %d\n[biWidth]: %d\n[biHeight]: %d\
\n[biPlanes]: %d\
\n[biBitCount]: %d\n[biCompression]: %d\
\n[biSizeImage]: %d\n[biXPelsPerMeter]: %d\
\n[biYPelsPerMeter]: %d\n[biClrUsed]: %d\
\n[biClrImportant]: %d\n",
bmp_i.biSize, bmp_i.biWidth, bmp_i.biHeight, bmp_i.biPlanes,
bmp_i.biBitCount, bmp_i.biCompression, bmp_i.biSizeImage,
bmp_i.biXPelsPerMeter, bmp_i.biYPelsPerMeter, bmp_i.biClrUsed, bmp_i.biClrImportant);
NUM_x = bmp_i.biWidth;
NUM_y = bmp_i.biHeight > 0 ? bmp_i.biHeight : -1 * bmp_i.biHeight;
// 分配buff
bmp_data = (unsigned char *)malloc(sizeof(char) * bmp_i.biSizeImage);
bmp_data_F = (struct _complex *)malloc( sizeof(struct _complex) * bmp_i.biSizeImage / 3);
memset(bmp_data, 0x00, bmp_i.biSizeImage);
fread(bmp_data, 1, bmp_i.biSizeImage, fp);
fclose(fp);
// 将 RGB 图转化为YUV灰度图， bmp 图像是小端存储，需要存 BGR
for(i = 0; i<NUM_y * NUM_x; i++){
b = bmp_data[3*i];
g = bmp_data[3*i+1];
r = bmp_data[3*i+2];
y = ( 66 * r + 129 * g + 25 * b + 128)/256 + 16 ;
//u = (unsigned char)( ( -38 * r - 74 * g + 112 * b + 128) >> 8) + 128 ;
//v = (unsigned char)( ( 112 * r - 94 * g - 18 * b + 128) >> 8) + 128 ;
if(y <= 0) y = 0;
else if(y >= 255) y = 255;
bmp_data[i] = (unsigned char)y;
bmp_data_F[i].x = (double)y;
bmp_data_F[i].y = 0;
}
// 如果该值 >0 说明bmp是倒存储的，需要对它做一个 Y镜像翻转
if(bmp_i.biHeight > 0){
for(i = 0; i<NUM_y; i++)
for(j = 0; j<NUM_x; j++){
bmp_data_F[i*NUM_x + j].x = (double)bmp_data[(NUM_y-i)*NUM_x + j];
}
}
// 将Y分量灰度图保存图片
//fp = fopen("orgin.yuv", "wb+");
//fwrite(bmp_data, 1, bmp_i.biSizeImage/3, fp);
//fclose(fp);
free(bmp_data);
}
void save_fft_yuv_file(char *url, struct _complex * buff)
{
int i, size = NUM_x * NUM_y;
bmp_data = (unsigned char *)malloc(size);
for(i = 0; i<size; i++){
bmp_data[i] = (unsigned char )sqrt(buff[i].x * buff[i].x + buff[i].y * buff[i].y);
}
FILE *fp = fopen(url, "wb+");
fwrite(bmp_data, 1, size , fp);
fclose(fp);
free(bmp_data);
}
int main(void)
{
open_bmp("少司命1024_512.bmp"); // 打开bmp图片，读取图片数据
fft(bmp_data_F, bmp_data_F, 1); // 对图片做傅里叶变换
save_fft_yuv_file("少司命1024_512_fft.yuv", bmp_data_F); // 将变换结果保存为 yuv 图片
fft(bmp_data_F, bmp_data_F, -1); // 对图片做傅里叶逆变换
save_fft_yuv_file("少司命1024_512_ifft.yuv", bmp_data_F); // 将逆变换结果保存为 yuv 图片
printf("\n\n程序运行结束\n");
return 0;
}
二、运行结果
2.1 《少司命1024_512.bmp》
文件《少司命1024_512.bmp》信息：
[bfType]: BM
[bfSize]: 1572918
[bfOffBits]:54
[biSize]: 40
[biWidth]: 1024
[biHeight]: 512
[biPlanes]: 1
[biBitCount]: 24
[biCompression]: 0
[biSizeImage]: 1572864
[biXPelsPerMeter]: 0
[biYPelsPerMeter]: 0
[biClrUsed]: 0
[biClrImportant]: 0
傅里叶变换结束，耗时0.854000s, NUM=1024 x 512
傅里叶逆变换结束，耗时0.872000s, NUM=1024 x 512
2.2 《test_少司命_256.bmp》
文件《test_少司命_256.bmp》信息：
[bfType]: BM
[bfSize]: 196662
[bfOffBits]:54
[biSize]: 40
[biWidth]: 256
[biHeight]: 256
[biPlanes]: 1
[biBitCount]: 24
[biCompression]: 0
[biSizeImage]: 196608
[biXPelsPerMeter]: 0
[biYPelsPerMeter]: 0
[biClrUsed]: 0
[biClrImportant]: 0
傅里叶变换结束，耗时0.078000s, NUM=256 x 256
傅里叶逆变换结束，耗时0.076000s, NUM=256 x 256
结果如下图：
分别为：变换前的BMP图、经过傅里叶变换后的频域图、经过傅里叶逆变换后的实域图。
2.3 《test_长方形.bmp》
2.4 《test_正方形.bmp》
2.5 《test_圆形.bmp》
2.6 《test_三角形.bmp》
2.7 《test_三角形&圆形.bmp》
可以看出，频率下各个频率是叠加的。
874
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
