fft_ifft.cpp

FFT和IFFT的轉換

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fft_ifft.h"

CFastFourier::CFastFourier()
{

}

CFastFourier::~CFastFourier()
{

}

void CFastFourier::swap(double* a, double* b)
{
    double t;
    
	 t	= *a;
    *a	= *b;
	*b	=  t;
}

void CFastFourier::bitrp(double* xreal, double* ximag, int n)
{
 	int i, j, a, b, p;

	for (i = 1, p = 0; i < n; i *= 2)
		p++;
    	
	for (i = 0; i < n; i ++)
    {
		a = i;
		b = 0;
		
		for (j = 0; j < p; j ++)
        {
			b = (b << 1) + (a & 1);	// b = b * 2 + a % 2;
			a >>= 1;				// a = a / 2;
        }
		
		if ( b > i)
        {
			swap(&xreal[i], &xreal[b]);
			swap(&ximag[i], &ximag[b]);
        }
     }
}


void CFastFourier::FFT(double* xreal, double* ximag, int n)
{
	double	wreal [N_ARRAY / 2], wimag [N_ARRAY / 2], treal, timag, ureal, uimag;
	double	arg;
    int		m, k, j, t, index1, index2;
     
    bitrp (xreal, ximag, n);
      
    arg			= -2 * PI / n;
    treal		= cos(arg);
    timag		= sin(arg);
    wreal [0]	= 1.0;
    wimag [0]	= 0.0;
     
	for (j = 1; j < n / 2; j ++)
    {
		wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
		wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
	}
 
    for (m = 2; m <= n; m *= 2)
    {
		for (k = 0; k < n; k += m)
        {
			for (j = 0; j < m / 2; j ++)
            {
				index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
    }
}

void CFastFourier::IFFT(double* xreal, double* ximag, int n)
{
	double	wreal [N_ARRAY / 2], wimag [N_ARRAY / 2], treal, timag, ureal, uimag;
	double	arg;
    int		m, k, j, t, index1, index2;
     
    bitrp (xreal, ximag, n);
    
    arg			= 2 * PI / n;
    treal		= cos(arg);
    timag		= sin(arg);
    wreal[0]	= 1.0;
    wimag[0]	= 0.0;
    
	for (j = 1; j < n / 2; j ++)
	{
		wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
        wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    }
 
    for (m = 2; m <= n; m *= 2)
	{
        for (k = 0; k < n; k += m)
		{
			for (j = 0; j < m / 2; j ++)
            {
				index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
	}
 
	for (j=0; j < n; j ++)
    {
		xreal[j] /= n;
        ximag[j] /= n;
	}
}

void CFastFourier::roundoff(double* zreal, int n)
{
	int i;
	
	for(i=0;i<n;i++)
		zreal[i] = (double)((int)(zreal[i] + 0.5));
}

void CFastFourier::mularray(double* xreal, double* ximag, double* yreal, double* yimag, double* zreal, double* zimag, int n)
{
	int i;
	
	for(i=0;i<n;i++)
	{
		zreal[i] = 	xreal[i] * yreal[i] - ximag[i] * yimag[i];
		zimag[i] =	xreal[i] * yimag[i] + yreal[i] * ximag[i];
	}
}

void CFastFourier::fft_conv(double* xreal, double* ximag, double* yreal, double* yimag, double* zreal, double* zimag, int n)
{
	FFT(xreal, ximag, n);
	FFT(yreal, yimag, n);

	mularray(xreal, ximag, yreal, yimag, zreal, zimag, n);
			
	IFFT(zreal, zimag, n);
}