fft.cpp

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

#include "dip.h"

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

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

void TDipWindow::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 TDipWindow::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 TDipWindow::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 TDipWindow::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 TDipWindow::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;
	}
}

void TDipWindow::tfft(int m,int n)
{
	int i,j;

	root(n);
	p = (COMPLEX *) new COMPLEX[n];
	f = (COMPLEX *) new COMPLEX[n];

	for ( i = 0 ; i < m ; i++ )
	{
		for ( j = 0 ; j < n ; j++ )
			p[j] = SpaceField[i][j];

		fft(n);

		for ( j = 0 ; j < n ; j++ )
			FreqField[i][j] = f[j];
	}

	delete p;
	delete f;

	for ( i = 0 ; i < m ; i++ )
		for ( j = 0 ; j < n ; j++ )
			SpaceField[i][j] = FreqField[i][j];

	root(m);
	p = (COMPLEX *) new COMPLEX[m];
	f = (COMPLEX *) new COMPLEX[m];

	for ( j = 0 ; j < n ; j++ )
	{
		for ( i = 0 ; i < m ; i++ )
			p[i] = SpaceField[i][j];

		fft(m);

		for ( i = 0 ; i < m ; i++ )
			FreqField[i][j] = f[i];
	}

	delete p;
	delete f;
}

void TDipWindow::ntfft(int m,int n)
{
	int i,j;

	root(n);
	p = (COMPLEX *) new COMPLEX[n];
	f = (COMPLEX *) new COMPLEX[n];

	for ( i = 0 ; i < m ; i++ )
	{
		for ( j = 0 ; j < n ; j++ )
			f[j] = FreqField[i][j];

		nfft(n);

		for ( j = 0 ; j < n ; j++ )
			SpaceField[i][j] = p[j];
	}

	delete p;
	delete f;

	for ( i = 0 ; i < m ; i++ )
		for ( j = 0 ; j < n ; j++ )
			FreqField[i][j] = SpaceField[i][j];

	root(m);
	p = (COMPLEX *) new COMPLEX[m];
	f = (COMPLEX *) new COMPLEX[m];

	for ( j = 0 ; j < n ; j++ )
	{
		for ( i = 0 ; i < m ; i++ )
			f[i] = FreqField[i][j];

		nfft(m);

		for ( i = 0 ; i < m ; i++ )
			SpaceField[i][j] = p[i];
	}

	delete p;
	delete f;
}