fft.c

超声波成像算法

#include <stdio.h>
#include <math.h> 
#include "Reconst.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846264338
#endif

void FFT(int nsamp,double *x, double *y, double f)	/* FFT */
	/* if f=-1.0 then fft */
	/* if f= 1.0 then rft */
{
//        char    ct;
	int mod=nsamp;
	int	i, i0, i1, j, l1, ns, k, arg;
	
	double	s, c, sc, x1, y1, t;

	#if 0
	   printf("fft start\n");
	   for(i=0;i<8;i++) printf("%d : %f\n",i,x[i]);
    #endif

	ns = mod/2;
	sc = 2*M_PI/(double)mod;

	while(ns >= 1) {
		arg = 0;
		for(l1 = 0; l1 < mod; l1 += (2*ns)) {
			k = mod/4;
			c = cos(sc*(double)arg);
			s = sin(f*sc*(double)arg);
			for(i0 = l1; i0 < l1+ns; i0++) {
				i1 = i0+ns;
				x1 = x[i1]*c-y[i1]*s;
				y1 = y[i1]*c+x[i1]*s;
				
				x[i1] = x[i0]-x1;
				y[i1] = y[i0]-y1;
				x[i0] = x[i0]+x1;
				y[i0] = y[i0]+y1;
			
			}
			while(k <= arg) {
				arg -= k;
				k /= 2;
				if(k ==0) break;
			}
			arg += k;
		}
		ns /= 2;
	}

	#if 0
	   printf("fft middle\n");
	   for(i=0;i<8;i++) printf("%d : %f\n",i,x[i]);
    #endif


	if(f < 0.0)
	for(i = 0; i < mod; i++) {
		x[i] /= (double)mod;
		y[i] /= (double)mod;
	}

	j = 0;
	for(i = 0; i < mod-1; i++) {
		if(i <= j) {
			t = x[i];	x[i] = x[j];	x[j] = t;
			t = y[i];	y[i] = y[j];	y[j] = t;
		}
		k = mod/2;
		while(k <= j) {
			j = j-k;
			k /= 2;
		}
		j = j+k;
	}

	#if 0
	   printf("fft end\n");
	   for(i=0;i<8;i++) printf("%d : %f\n",i,x[i]);
    #endif

}