fft_flt.h

用systemC制作的快速付立叶变换

// fft_flt.h
// Pierre Bomel - LESTER/UBS - janvier 2005
//    Formation chez France Telecom R&D
#include <systemc.h>
#include <math.h>

template <int n, int log2n, bool inv, bool trace_me>
// n = number of points
// log2n = log2(n)
// inverse = inverse FFT ON/OFF
// trace_me = trace ON/OFF
SC_MODULE (fft_flt){
	// channels
	sc_in<double> sreals[n], simags[n];
	sc_inout<double> freals[n], fimags[n];
	// trace
	ofstream coutr;
	ofstream couti;

	// data to process
	double coeffreals[n], coeffreal, treals[n], topreal, botreal; // real parts
	double coeffimags[n], coeffimag, timags[n], topimag, botimag; // imaginary parts
	// FFT indexes
	int Kn, Bp, Np, P, Np2, BaseT, b, BaseB, k;

	void init_coeffs(){
		double sign = inv ? 1.0 : -1.0;
		double delta = (M_PI * 2.0) / double(n);
		double angle = 0;
		for(int i = 0; i < n; i++){
			coeffreals[i] = cos(angle);
			//cout << i << " cos " << coeffreals[i] << endl;
			coeffimags[i] = sign * sin(angle);
			//cout << i << " sin " << coeffimags[i] << endl;
			angle += delta;
		}
		Kn = n;
	}
	double Tr(int np, int i) {return coeffreals[i * Kn/np];}
	double Ti(int np, int i) {return coeffimags[i * Kn/np];}

	int bit_reverse(int index)
	{
		int inverse = 0;
		for(int i = 0; i < log2n; i++){
			inverse <<= 1;
			inverse |= index & 0x1;
			index >>= 1;
		}
		return inverse;
	}

	void fft(){

		// copy input channels values into array in bit reverse order
		for(int i = 0; i < n; i++){
			int br = bit_reverse(i);
			treals[br] = sreals[i].read();
			timags[br] = simags[i].read();
			//cout << i << "(" << br << ") " << treals[br] << endl;
			//cout << i << "(" << br << ") " << timags[br] << endl;
		}

		// run DIT radix2 FFT
		Bp = 1 << (log2n-1);
		Np = 2;
		for(P = 0; P < log2n; P++){
			Np2 = Np >> 1;
			BaseT = 0;
			for(b = 0; b < Bp; b++){
				BaseB = BaseT + Np2;
				for(k = 0; k < Np2; k++){
					int l1 = BaseT+k;
					int l2 = BaseB+k;
					topreal = treals[l1];
					//cout << "topreal " << topreal << endl;
					topimag = timags[l1];
					//cout << "topimag " << topimag << endl;
					coeffreal = Tr(Np, k);
					//cout << "coeffreal " << coeffreal << endl;
					coeffimag = Ti(Np, k);
					//cout << "coeffimag " << coeffimag << endl;
					botreal = treals[l2]*coeffreal - timags[l2]*coeffimag;
					//cout << "botreal " << botreal << endl;
					botimag = treals[l2]*coeffimag + timags[l2]*coeffreal;
					//cout << "botimag " << botimag << endl;
					treals[l1] = topreal + botreal;
					//cout << l1 << " " << treals[l1] << endl;
					timags[l1] = topimag + botimag;
					//cout << l1 << " " << timags[l1] << endl;
					treals[l2] = topreal - botreal;
					//cout << l2 << " " << treals[l2] << endl;
					timags[l2] = topimag - botimag;
					//cout << l2 << " " << timags[l2] << endl;
				}
				BaseT = BaseT + Np;
			}
			Bp = Bp >> 1;
			Np = Np << 1;
		}
		
		// scaling only for reverse FFT (like MATLAB does)
		if (inv){
			for(int i = 0; i < n; i++){
				treals[i] /= n;
				timags[i] /= n;
			}
		}

		// trace
		if (trace_me) {
			for(int i = 0; i < n; i++){
				coutr << i << " " << treals[i] << endl;
				couti << i << " " << timags[i] << endl;
			}
		}

		// copy array into output channels values
		for(int i = 0; i < n; i++){
			freals[i].write(treals[i]);
			fimags[i].write(timags[i]);
		}
	}

	SC_CTOR(fft_flt){
		if (trace_me) {
			char file_name[100];
			sprintf(file_name, "%sfft_flt_real", inv ? "i" : "");
			coutr.open(file_name, ios::out);
			sprintf(file_name, "%sfft_flt_imag", inv ? "i" : "");
			couti.open(file_name, ios::out);
		}
		init_coeffs();
		SC_METHOD(fft);
		for(int i = 0; i < n; i++){
			sensitive << sreals[i] << simags[i];
		}
		dont_initialize();
	}
};