testfft.cpp

很经典的数字图像处理讲义

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "temp.h"

COMPLEX *omega,*p,*f,*temp;

void add(COMPLEX *a,COMPLEX *b,COMPLEX *ret)
{
	ret->re = a->re + b->re;
	ret->im = a->im + b->im;
}

void sub(COMPLEX *a,COMPLEX *b,COMPLEX *ret)
{
	ret->re = a->re - b->re;
	ret->im = a->im - b->im;
}

void mul(COMPLEX *m,COMPLEX *n,COMPLEX *ret)
{
	float a,b,c;

	a = ( m->re - m->im ) * n->im;
	b = m->re * ( n->re - n->im );
	c = m->im * ( n->re + n->im );
	ret->re = a + b;
	ret->im = a + c;
}

int reverse(int t,int k)
{
	int i,j,m,s;

	m = t;
	s = 0;
	for ( i = 0 ; i < k ; i++ )
	{
		j = m / 2;
		s = s * 2 + ( m - j * 2 );
		m = j;
	}

	return s;
}

void root(int n)
{
	int j;

	omega->re = 1;
	omega->im = 0;
	omega[1].re = (float)cos(2 * PI / n);
	omega[1].im = -(float)sin(2 * PI / n);

	for ( j = 2 ; j < n/2 ; j++ )
		mul(&omega[1],&omega[j-1],&omega[j]);

	for ( j = n/2 ; j < n ; j++ )
	{
		omega[j].re = -omega[j-n/2].re;
		omega[j].im = -omega[j-n/2].im;
	}
}

void fft(int n)
{
	int s,k,m,l,nv,t,j;
	COMPLEX podd,ret;

	k = (int)(log(n) / log(2) + 0.5);
	nv = n;
	m = 1;
	for ( l = k-1 ; l >= 0 ; l-- )
	{
		for ( t = 0 ; t < m * nv ; t+=nv )
			for ( j = 0 ; j < nv/2 ; j++ )
			{
				s = (t+j) / (int)pow(2,l);
				s = reverse(s,k);

				ret = omega[s];
				mul(&ret,&p[t+j+nv/2],&podd);
				sub(&p[t+j],&podd,&p[t+j+nv/2]);
				add(&p[t+j],&podd,&p[t+j]);
			}
		m *= 2;
		nv /= 2;
	}

	for ( t = 0 ; t < n ; t++ )
	{
		s = reverse(t,k);
		f[t] = p[s];
	}
}

void nfft(int n)
{
	int i;

	for ( i = 0 ; i < n ; i++ )
		p[i] = f[i];

	fft(n);

	p[0] = f[0];
	for ( i = 1 ; i < n ; i++ )
		p[i] = f[n-i];
	for ( i = 0 ; i < n ; i++ )
	{
		p[i].re /= n;
		p[i].im /= n;
	}
}

int main()
{
	omega = new COMPLEX[8];
	p = new COMPLEX[8];
	f = new COMPLEX[8];
	temp = new COMPLEX[8];

	root(8);

	p[0].re = 1;
	p[1].re = 1;
	p[2].re = 2;
	p[3].re = 3;
	p[4].re = 0;
	p[5].re = 0;
	p[6].re = 0;
	p[7].re = 0;
	for ( int i = 0 ; i < 8 ; i++ )
		p[i].im = 0;

	fft(8);
	for ( i = 0 ; i < 8 ; i++ )
	{
		printf("%f  ,  %f\n",f[i].re,f[i].im);
		temp[i].re = f[i].re;
		temp[i].im = f[i].im;
	}
	printf("\n");

	p[0].re = 2;
	p[1].re = 3;
	p[2].re = 2;
	p[3].re = 4;
	p[4].re = 0;
	p[5].re = 0;
	p[6].re = 0;
	p[7].re = 0;
	for ( i = 0 ; i < 8 ; i++ )
		p[i].im = 0;

	fft(8);
	for ( i = 0 ; i < 8 ; i++ )
	{
		printf("%f  ,  %f\n",f[i].re,f[i].im);
		mul(&f[i],&temp[i],&f[i]);
	}
	printf("\n");

	for ( i = 0 ; i < 8 ; i++ )
	{
		printf("%f  ,  %f\n",f[i].re,f[i].im);
	}
	printf("\n");

	nfft(8);
	for ( i = 0 ; i < 8 ; i++ )
	{
		printf("%f  ,  %f\n",p[i].re,p[i].im);
	}
	printf("\n");

	return 0;
}