数学计算

开发平台：
C++ Builder

fft.cpp：源码内容
							#include "stdafx.h"
#include "fft.h"
#include "Complex.h"
#include<math.h>
#define M_PI 3.141592654		// The mathematical constant for PI
void FFT1D(Complex<double> *fx, Complex<double> *Fu, int twoK, int stride) {
	// This function performs the 1D fast fourier transform from input data fx,
	// it outputs the frequency coefficients in the array Fu.
	// Input:	fx - pointer to an array of data in complex number format
	//			twoK - the total number of data in fx
	//			stride - the number of data to be skipped when reading the next data in fx
	//					e.g. fx = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15], twoK = 4
	//						if stride = 1, the data to be processed is [0 1 2 3]
	//						if stride = 2, the data to be processed is [0 2 4 8]
	//						if stride = 4, the data to be processed is [0 4 8 12]
	// Output: Fu - pointer to an array buffer for transformed cofficients in complex number format
	// Code to be added here ...
	Complex<double> *F_even,*F_odd, e;
	int u,K;
	if(twoK==2)
	{
		Fu[0]=(fx[0]+fx[1*stride])*0.5;
		Fu[1]=(fx[0]-fx[1*stride])*0.5;
		return;
	}
	K=twoK/2;
	F_even=new Complex<double>[K];
	F_odd=new Complex<double>[K];
	FFT1D(fx,F_even,K,2*stride);
	FFT1D(fx+stride, F_odd, K, 2*stride);
	for(u=0;u<K;u++)
	{
		e.Exp(-M_PI*((double)u/(double)K));
		Fu[u]=(F_even[u]+e*F_odd[u])*0.5;
		Fu[u+K]=(F_even[u]-e*F_odd[u])*0.5;
	}
	delete []F_even;
	delete []F_odd;
}
void FFT2D(Complex<double> *f, Complex<double> *F, int width, int height) {
	// This function performs the 2D Fast Fourier Transform for a given image f
	// You can safely assume that the width = 2^m, height = 2^n for some integers m and n,
	//		and the input image are already properly padded
	// Input: f - An 2D Image represented in Complex Numbers
	// Output: F - The transformed coefficients, also represented in Complex Numbers
	// Code to be added here ...
	int i,j;
	Complex<double> **s = new Complex<double>*[height];
	for(int j=0; j<height; j++){
	   s[j] =new Complex<double>[width];
	}	
	for(j=0;j<height;j++)
		for(i=0;i<width;i++)
			s[j][i]=f[(j*width)+i];
	Complex<double> *row;
	row=new Complex<double>[width];
	for(j=0;j<height;j++)
	{
		for(i=0;i<width;i++)
			row[i]=s[j][i];
		FFT1D(row,row,width,1);
		for(i=0;i<width;i++)
			s[j][i]=row[i];
	}
	Complex<double> *col;
	col=new Complex<double>[height];
	for(i=0;i<width;i++)
	{
		for(j=0;j<height;j++)
			col[j]=s[j][i];
		FFT1D(col,col,height,1);
		for(j=0;j<height;j++)
			s[j][i]=col[j];
	}
	for(j=0;j<height;j++)
		for(i=0;i<width;i++)
			F[(j*width)+i]=s[j][i];
	delete []row;
	delete []col;
	delete s;
	
}
void IFFT2D(Complex<double> *F, Complex<double> *f, int width, int height) {
	// This function performs the Inverse 2D Fast Fourier Transform for a given image f (in one color channel R, G or B)
	// You can safely assume that the width = 2^m, height = 2^n for some integers m and n
	// Input: F - The transformed coefficients, also represented in Complex Numbers
	// Output: f - An 2D Image represented in Complex Numbers
	// Code to be added here ...
	int i;
	for(i=0;i<width*height;i++)
		{
			F[i]=Complex<double>(F[i].Real(),-F[i].Imag());
		}
	FFT2D(F,f,width,height);
	for(i=0;i<width*height;i++)
		{
		f[i]=Complex<double>(f[i].Real(),-f[i].Imag());
		f[i]=f[i]*width*height;
		}
}
void ScaleFFT(Complex<double> *F, Complex<double> *PS, int width, int height) {
	// This function takes a 2D coefficients in F, and compute the power spectrum of F.
	// This function also scales the power spectrum, such that its values lie within the range [0, 255]
	//		for proper display.
	// Input:	F - 2D coefficients of the FFT image
	//			width - width of the F
	//			height - height of the F
	// Output:	PS - scaled power spectrum (with value normalized to the range of [0, 255])
	// Code to be added here ...
	
int i;
double min,max,diff;
for(i=0;i<width*height;i++)
	{
		PS[i] = Complex <double>(pow((F[i].Mod()*F[i].Mod()),0.2));
	}
min=max=PS[0].Real();
for(i=1;i<width*height;i++)
{
if(min > PS[i].Real())
	min = PS[i].Real();
if(max < PS[i].Real())
	max =PS[i].Real();
}
diff = max-min;
for(i=0;i<width*height;i++)
	{
	PS[i] = Complex <double>(((PS[i].Real()-min)/diff)*255.0);
	}
}
void BLPF(Complex<double> *F, Complex<double> *LF, int width, int height) {
	// This function takes a 2D coefficients in F, and perform a second other ButterWorth
	//		low pass filtering on the coefficents. The filtered coefficients are put back into LF
	// Input:	F - 2D coefficients of the FFT image
	//			width - width of the F
	//			height - height of the F
	// Output:	LF - ButterWorth filtered coefficients
	// The coded below is a low pass filter implementation, however, it cannot filter the image
	// properly. Modify the code below such that the filter becomes a butterworth filter of order 2.
	/*int i, j ;
	int Df ;
	int centerX, centerY ;
	centerX = width / 2 ;
	centerY = height / 2 ;
	Df = 32 ;
	for(j = 0; j < height; j++) {
		for(i = 0; i < width; i++) {
			if((i - centerX) * (i - centerX) + (j - centerY) * (j - centerY) > Df * Df) {
				LF[i + j * width] = 0.0 ;
			} else {
				LF[i + j * width] = F[i + j * width] ;
			}
		}
	}*/
	
	int i, j,D;
	double T,H;
	int Df ;
	int centerX, centerY ;
	
	centerX = width / 2 ;
	centerY = height / 2 ;
	Df = 32 ;
	for(j = 0; j < height; j++) {
		for(i = 0; i < width; i++) {
			D=(i - centerX) * (i - centerX) + (j - centerY) * (j - centerY);
			T=(D/(Df*Df))*(D/(Df*Df));
			H=1/(1+T);
			LF[i+j*width]=F[i+j*width]*H;
			}
	}
}

						