2-FFT变换

最新推荐文章于 2021-09-13 03:11:30 发布

CUC_ENGINEER

最新推荐文章于 2021-09-13 03:11:30 发布

阅读量1.1k

点赞数 1

分类专栏：数字图像处理

本文链接：https://blog.csdn.net/CUC_ENGINEER/article/details/105598656

版权

数字图像处理专栏收录该内容

2 篇文章 0 订阅

订阅专栏

实验原理

2-D傅里叶变换

根据2-D傅里叶变换的可分离性，可以将二维的傅里叶变换分解为两次一维的傅里叶变换实现，一维的傅里叶变换可以通过快速傅里叶算法实验，从而大大减少计算量。

##2-D傅里叶逆变换

二维傅里叶变换的逆变换可以利用二维傅里叶正变换来实现。

正变换：
$F（u，v)=\frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)e^{-j2\pi(ux/M+vy/N)}$

反变换:
$\begin{aligned} F（x，y)&=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1}f(u,v)e^{j2\pi(ux/M+vy/N)} \\&=MN*\left(\frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)^*e^{-j2\pi(ux/M+vy/N)}\right)^* \end{aligned}$

实验步骤

2-Dfft

预处理：由于使用了FFT算法，所以必须保证输入的长和宽为2的整次幂，不够补0
先通过输入图像的每一行计算一维变换
再沿中间结果的每一列计算一维变换

注意：由于读入数据使按行排列，所以计算列向量时要进行转置
再经过转置得到处理后的结果

2-Difft

二维傅里叶逆变换相对简单

取共轭
2-Dfft变换
取共轭，乘系数
去除FFT时补入的零点

实验结果

可以看出两者并没有太大区别：
原图像

原图像

处理后的图像

实现代码

/*************************************主函数***************************************/
#include <stdio.h> 
#include <stdlib.h>
#include <malloc.h>
#include <math.h>   
#include <assert.h>
#include "fft.h"
#include "fft_2d.h"
#define WIDTH   464
#define HEIGHT 538
#define eps   1e-10

int main(int argc, char *argv[])
{
	unsigned char *yBuf;
	yBuf=(unsigned char *)malloc(WIDTH*HEIGHT);
	FILE* AFile =NULL;
	FILE* BFile =NULL;
	AFile=fopen(argv[1],"rb");
	BFile=fopen(argv[2],"wb");
	if(AFile==NULL)
	{
		printf("cannot find AFile \n");
		exit(1);
	}
	if(BFile==NULL)
	{
		printf("cannot find BFile \n");
		exit(1);
	}
	if (!feof(AFile))
	{
		if(!fread(yBuf,1,WIDTH*HEIGHT,AFile))
		{
			printf("cannot read y%d ");
		}
	}
	int width=mod2(WIDTH);
	int height=mod2(HEIGHT);
	double *A_re, *A_im; 
	int N=width*height;
	A_re = (double*)malloc(sizeof(double)*N); 
	A_im = (double*)malloc(sizeof(double)*N);
	assert(A_re != NULL && A_im != NULL); 
	for (int n=0; n<N; n++) 
	{
		A_re[n] = 0;
		A_im[n] = 0;
	}
	//将图片中读出的char转成double
	for (int n=0;n<N;n++ )
		if (n<WIDTH*HEIGHT)
				A_re[n]=(double)yBuf[n];
	fft_2d(WIDTH,HEIGHT,A_re,A_im);
	ifft_2d(WIDTH,HEIGHT,A_re,A_im);
	//将变换后的的double转成char
	for (int n=0;n<WIDTH*HEIGHT;n++ )
	{
		if (A_re[n]>=256)
		{
			printf("转换溢出");
			getchar();
		}
		else
			yBuf[n]=(unsigned char)A_re[n];
	}
	fwrite(yBuf,1,WIDTH*HEIGHT,BFile);
	for (int n=0;n<WIDTH*HEIGHT;n++ )
		yBuf[n]=128;//填充U，V
	fwrite(yBuf,1,WIDTH*HEIGHT,BFile);
	fwrite(yBuf,1,WIDTH*HEIGHT,BFile);
	free(A_im);
	A_im=NULL;
	free(A_re);
	A_re=NULL;
	free(yBuf);
	yBuf=NULL;
	fclose(AFile);
	fclose(BFile);
	return 0;
}

/***********************************二维FFT变换*************************************/
#include <stdio.h> 
#include <stdlib.h>
#include <malloc.h>
#include <math.h>   
#include <assert.h>
#include "fft.h"
#include "fft_2d.h"

void fft_2d(int W,int H,double *A_re, double *A_im)
{
	addzero_to2( W , H,A_re, A_im);//补零至2的整次幂
	int width=mod2(W);//向上取整到2的整次幂
	int height=mod2(H);
	int N=width*height;
	for (int n=0;n<height;n++)
		fft_1d(width, A_re+n*width, A_im+n*width); //水平方向做一维FFT变换
	double *temp_re,*temp_im;
	temp_re = (double*)malloc(sizeof(double)*N);
	temp_im = (double*)malloc(sizeof(double)*N);
	/** 矩阵转置注意不要越界**/
	for (int n=0;n<width;n++)
		for (int m=0;m<height;m++)
		{
			temp_re[n*height+m]=A_re[m*width+n];
			temp_im[n*height+m]=A_im[m*width+n];
		}
	for (int m=0;m<width;m++)
		fft_1d(height,temp_re+m*height, temp_im+m*height); //垂直方向做一维FFT变换
	for (int m=0;m<height;m++)
	{
		for (int n=0;n<width;n++)
		{
			A_re[m*width+n]=temp_re[n*height+m];
			A_im[m*width+n]=temp_im[n*height+m];
		}
	}
	free(temp_re);
	temp_re=NULL;
	free(temp_im);
	temp_im=NULL;
}

void ifft_2d(int W,int H,double *A_re, double *A_im)//二维iFFT变换
{
	int width=mod2(W);//向上取整到2的整次幂
	int height=mod2(H);
	for (int n=0;n<height;n++)//取共轭
		for (int m=0;m<width;m++)
			A_im[n*width+m]=-A_im[n*width+m];
	fft_2d(width,height,A_re,A_im);//二维FFT变换
	for (int n=0;n<height;n++)//取共轭
		for (int m=0;m<width;m++)
		{
			A_im[n*width+m]=-A_im[n*width+m]/(width*height);
			A_re[n*width+m]=A_re[n*width+m]/(width*height);
		}
	minuszero_from2(W ,H,A_re , A_im );//去除尾端的0
}

void addzero_to2(int W ,int H,double *A_re, double *A_im)//补零至2的整次幂
{
	int width=mod2(W);
	int height=mod2(H);
	int N=width*height;
	double *temp_re,*temp_im;
	temp_re = (double*)malloc(sizeof(double)*N);
	temp_im = (double*)malloc(sizeof(double)*N);
	for (int n=0; n<N; n++) 
	{
		temp_re[n] = 0;
		temp_im[n] = 0;
	}
	for (int i=0; i<H;i++)
	{
		for (int j=0;j<W;j++)
		{
			temp_re[i*width+j]=A_re[i*W+j];
			temp_im[i*width+j]=A_im[i*W+j]; 
		}
	}
	for (int n=0; n<N; n++) 
	{
		A_re[n] = temp_re[n];
		A_im[n] = temp_im[n];
	}
	free(temp_re);
	temp_re=NULL;
	free(temp_im);
	temp_im=NULL;
}
void minuszero_from2(int W ,int H,double *A_re, double *A_im)//去除尾端的0
{
	int width=mod2(W);
	int height=mod2(H);
	int N=W*H;
	double *temp_re,*temp_im;
	temp_re = (double*)malloc(sizeof(double)*N);
	temp_im = (double*)malloc(sizeof(double)*N);
	for (int n=0; n<N; n++) 
	{
		temp_re[n] = 0;
		temp_im[n] = 0;
	}
	for (int i=0; i<H;i++)
	{
		for (int j=0;j<W;j++)
		{
			temp_re[i*W+j]=A_re[i*width+j];
			temp_im[i*W+j]=A_im[i*width+j]; 
		}
	}
	for (int n=0; n<N; n++) 
	{
		A_re[n] = temp_re[n];
		A_im[n] = temp_im[n];
	}
	free(temp_re);
	temp_re=NULL;
	free(temp_im);
	temp_im=NULL;
}
int mod2(int a)//向上取整到2的整次幂
{
	for (a; (a & (a - 1)) ; a++);
	return a;
}

/***********************************一维FFT变换*************************************/
#include <stdio.h> 
#include <stdio.h> 
#include <stdlib.h>
#include <math.h>   
#include <assert.h>

#include "fft.h"

#define  M_PI 3.1415926


/***********************************************************************************
对N点序列进行fft变换：
1. A_re为该序列的实部，A_im为该序列的虚部；
2. 运算结果仍然存放在A_re和A_im 
3. 输入输出都是自然顺序

Note：
N必须为2的整数次幂
**********************************************************************************/
void fft_1d(int N , double *A_re, double *A_im) 
{ 
	double *W_re, *W_im;

	W_re = (double*)malloc(sizeof(double)*N/2); 
	W_im = (double*)malloc(sizeof(double)*N/2); 
	assert(W_re != NULL && W_im != NULL); 

	compute_W(N, W_re, W_im);     
	fft_butterfly(N, A_re, A_im, W_re, W_im);
	permute_bitrev(N, A_re, A_im); 

	free(W_re); 
	free(W_im);

}



/* W will contain roots of unity so that W[bitrev(i,log2n-1)] = e^(2*pi*i/n)
* n should be a power of 2
* Note: W is bit-reversal permuted because fft(..) goes faster if this is done.
*       see that function for more details on why we treat 'i' as a (log2n-1) bit number.
*/
void compute_W(int n, double *W_re, double *W_im)
{
	int i, br;
	int log2n = log_2(n);

	for (i=0; i<(n/2); i++)
	{
		br = bitrev(i,log2n-1); 
		W_re[br] = cos(((double)i*2.0*M_PI)/((double)n));  
		W_im[br] = -sin(((double)i*2.0*M_PI)/((double)n));  
	}

}


/* permutes the array using a bit-reversal permutation */ 
void permute_bitrev(int n, double *A_re, double *A_im) 
{ 
	int i, bri, log2n;
	double t_re, t_im;

	log2n = log_2(n); 

	for (i=0; i<n; i++)
	{
		bri = bitrev(i, log2n);

		/* skip already swapped elements */
		if (bri <= i) continue;

		t_re = A_re[i];
		t_im = A_im[i];
		A_re[i]= A_re[bri];
		A_im[i]= A_im[bri];
		A_re[bri]= t_re;
		A_im[bri]= t_im;
	}  
} 


/* treats inp as a numbits number and bitreverses it. 
* inp < 2^(numbits) for meaningful bit-reversal
*/ 
int bitrev(int inp, int numbits)
{
	int i, rev=0;
	for (i=0; i < numbits; i++)
	{
		rev = (rev << 1) | (inp & 1);
		inp >>= 1;
	}
	return rev;
}


/* returns log n (to the base 2), if n is positive and power of 2 */ 
int log_2(int n) 
{
	int res; 
	for (res=0; n >= 2; res++) 
		n = n >> 1; 
	return res; 
}

/* fft on a set of n points given by A_re and A_im. Bit-reversal permuted roots-of-unity lookup table
* is given by W_re and W_im. More specifically,  W is the array of first n/2 nth roots of unity stored
* in a permuted bitreversal order.
*
* FFT - Decimation In Time FFT with input array in correct order and output array in bit-reversed order.
*
* REQ: n should be a power of 2 to work. 
*
* Note: - See www.cs.berkeley.edu/~randit for her thesis on VIRAM FFTs and other details about VHALF section of the algo
*         (thesis link - http://www.cs.berkeley.edu/~randit/papers/csd-00-1106.pdf)
*       - See the foll. CS267 website for details of the Decimation In Time FFT implemented here.
*         (www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html)
*       - Also, look "Cormen Leicester Rivest [CLR] - Introduction to Algorithms" book for another variant of Iterative-FFT
*/
void fft_butterfly(int n, double *A_re, double *A_im, double *W_re, double *W_im) 
{
	double w_re, w_im, u_re, u_im, t_re, t_im;
	int m, g, b;
	int mt, k;

	/* for each stage */  
	for (m=n; m>=2; m=m>>1) 
	{
		/* m = n/2^s; mt = m/2; */
		mt = m >> 1;

		/* for each group of butterfly */ 
		for (g=0,k=0; g<n; g+=m,k++) 
		{
			/* each butterfly group uses only one root of unity. actually, it is the bitrev of this group's number k.
			* BUT 'bitrev' it as a log2n-1 bit number because we are using a lookup array of nth root of unity and
			* using cancellation lemma to scale nth root to n/2, n/4,... th root.
			*
			* It turns out like the foll.
			*   w.re = W[bitrev(k, log2n-1)].re;
			*   w.im = W[bitrev(k, log2n-1)].im;
			* Still, we just use k, because the lookup array itself is bit-reversal permuted. 
			*/
			w_re = W_re[k];
			w_im = W_im[k];

			/* for each butterfly */ 
			for (b=g; b<(g+mt); b++) 
			{
				/* t = w * A[b+mt] */
				t_re = w_re * A_re[b+mt] - w_im * A_im[b+mt];
				t_im = w_re * A_im[b+mt] + w_im * A_re[b+mt];

				/* u = A[b]; in[b] = u + t; in[b+mt] = u - t; */
				u_re = A_re[b];
				u_im = A_im[b];
				A_re[b] = u_re + t_re;
				A_im[b] = u_im + t_im;
				A_re[b+mt] = u_re - t_re;
				A_im[b+mt] = u_im - t_im;

			}
		}
	}
}

CUC_ENGINEER

关注

1
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
2-FFT变换

实验原理2-D傅里叶变换根据2-D傅里叶变换的可分离性，可以将二维的傅里叶变换分解为两次一维的傅里叶变换实现，一维的傅里叶变换可以通过快速傅里叶算法实验，从而大大减少计算量。##2-D傅里叶逆变换二维傅里叶变换的逆变换可以利用二维傅里叶正变换来实现。正变换：F（u，v)=1MN∑x=0M−1∑y=0N−1f(x,y)e−j2π(ux/M+vy/N)F（u，v)=\frac{1}{MN...
复制链接

扫一扫