tktfft.cpp

用Matlab与C的混合编程写的FFT（快速傅立叶变换）算法

// Function calculate fft of an array x
// X=tktfft(x,N),if you want to specify N
// or X=tktfft(x),when using default N=512
// Ather , 3/27

#include<math.h>
#include "mex.h"
#define pi 3.141592653


void mexFunction( int nlhs,mxArray *plhs[],int nrhs,mxArray *prhs[]){
    void cal(mxArray *,int,mxArray *);
	void evenHalf(mxArray *,int,mxArray *);
	void oddHalf(mxArray *,int,mxArray *);
// Error checking
	if( nrhs > 2 )
		mexErrMsgTxt("Too many input arguments!");
    else if( nlhs >1 )
		mexErrMsgTxt("Too many output arguments!");
	else if( nlhs<1 )
		mexErrMsgTxt("Too few output arguments!");
	else if( !mxIsNumeric(prhs[0]) )
			mexErrMsgTxt("Input array should be numeric!");
	else if( mxGetM(prhs[0]) != 1 )
		mexErrMsgTxt("Sorry,I can just handle row vectors.. (:");

    mxArray * input=prhs[0];
    int incol=mxGetN(input);		// how many collums we get
	int N;
	if( nrhs==1)		//only one input parameter
		N=incol;
	else 
		N=(int) mxGetScalar(prhs[1]);			// parameters start form 0

	plhs[0]=mxCreateDoubleMatrix(1,incol,mxCOMPLEX);
    cal(input,N,plhs[0]);
    return;
}

void evenHalf(mxArray * input,int incol,mxArray * out){
	int k;
    double * rout = mxGetPr(out);   
    double * rin = mxGetPr(input);
    double * iout,* iin;
    for(k=0;k<incol/2;k++)
        rout[k] = rin[k*2];
    if(mxIsComplex(input) ){
        iout = mxGetPi(out);
         iin = mxGetPi(input);
        for(k=0;k<incol/2;k++)
            iout[k] = iin[k*2];
    }
    return;
}

void oddHalf(mxArray * input,int incol,mxArray * out){
	int k;
    double * rout = mxGetPr(out);  
    double * rin = mxGetPr(input);
    double * iout,* iin;
    for(k=0;k<incol/2;k++)
       rout[k] = rin[k*2+1];
    if(mxIsComplex(input) ){
        iout = mxGetPi(out);
        iin = mxGetPi(input);
        for(k=0;k<incol/2;k++)
            iout[k] = iin[k*2+1];
    }
    return;
}

void cal(mxArray *input,int N,mxArray *out){
//	calculate subroutine
	int k;
	int incol=N;		// how many collums we get
	if(incol==1){                   // the exit of iteration
//		out=mxDuplicateArray(input);
		double * rout = mxGetPr(out);
		double * rinput = mxGetPr(input);
		rout[0] = rinput[0];
		if( mxIsComplex(input) ){
			double * iout = mxGetPi(out);
			double * iinput = mxGetPi(input);
			iout[0] = iinput[0];
		}
		return;
	}
	else{
        mxArray *part1 = mxCreateDoubleMatrix(1,incol/2,mxCOMPLEX);  evenHalf(input,incol,part1);
        mxArray *part2 = mxCreateDoubleMatrix(1,incol/2,mxCOMPLEX);  oddHalf(input,incol,part2);

        mxArray *FFT_part1 = mxCreateDoubleMatrix(1,incol/2,mxCOMPLEX);
        mxArray *FFT_part2 = mxCreateDoubleMatrix(1,incol/2,mxCOMPLEX);
// call this program recursively        
		cal(part1,N/2,FFT_part1);
		cal(part2,N/2,FFT_part2);
// get real part and imaginary part of array        
        double *rFFT_part1 = mxGetPr(FFT_part1);
        double *iFFT_part1 = mxGetPi(FFT_part1);
        double *rFFT_part2 = mxGetPr(FFT_part2); 
        double *iFFT_part2 = mxGetPi(FFT_part2);
        double *rplhs = mxGetPr(out);
        double *iplhs = mxGetPi(out);

		for(k=0;k<incol/2;k++){
			rplhs[k]=rFFT_part1[k]+rFFT_part2[k]*cos(2*pi/incol*k) + iFFT_part2[k]*sin(2*pi/incol*k);
            iplhs[k]=iFFT_part1[k]+iFFT_part2[k]*cos(2*pi/incol*k) - rFFT_part2[k]*sin(2*pi/incol*k);
            rplhs[k+incol/2]=rFFT_part1[k]-rFFT_part2[k]*cos(2*pi/incol*k) - iFFT_part2[k]*sin(2*pi/incol*k);
            iplhs[k+incol/2]=iFFT_part1[k]-iFFT_part2[k]*cos(2*pi/incol*k) + rFFT_part2[k]*sin(2*pi/incol*k);
		}// end for
		return;
	}// end else
}