fft_2d.c

数字信号处理实验

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define PI        3.1415926535

/********$B%W%m%H%?%$%W@k8@(B********/
void FFT1(double *re_part, double *im_part,
	  int num_of_data, int flag);

void FFT2(double *redata, double *imdata,
	  int width, int height, int flag)
/*

double$B7?$KJQ49$7$?#2<!85%G!<%?$N%U!<%j%(JQ49$r9T$&(B
$B2hA|$N>l9g$O;vA0$K(Bdouble$B7?$KJQ49$9$kI,MW$"$j(B
flag$B$r(B 1$B$G<B9T$9$k$H%U!<%j%(JQ49(B
flag$B$r(B-1$B$G<B9T$9$k$H%U!<%j%(5UJQ49(B

$BI}(B,$B9b$5$r<u$1$k!#(B

$B;29M(B:
$B#C8@8l$K$h$k2hA|=hM}F~Lg(B $B><98F2(B $B0B5o1!(B $BLT(B $BCx(B
IDBN4-7856-3124-4 3000yen

$B0z?t$r2~B$:Q$_(B

*/

{
  int i,j,mat_size;
  double *re, *im;

  mat_size = width * height;

  re = malloc(width*sizeof(double));
  im = malloc(width*sizeof(double));

  /*  
  printf("%d, %d, %d\n",width, height, flag);
  */

  if(flag == 1){
    printf("FFT....\n");
  }
  else if(flag == -1){
    printf("Inv. FFT....\n");
  }
  else   {
    exit(0);
  }

  for(i=0;i<height;i++)
    {
      for(j=0;j<width;j++)
	{
          re[j]=redata[i*width + j];
	  im[j]=imdata[i*width + j];
	}
      
      FFT1(re, im, width, flag);
      
      for(j=0;j<width;j++)
	{
          redata[i*width + j]= re[j];
	  imdata[i*width + j]= im[j];
	}
    }

  for(i=0;i<width;i++)
    {
      for(j=0;j<height;j++)
	{
          re[j]=redata[j*width + i];
	  im[j]=imdata[j*width + i];
	}
      
      FFT1(re, im, height, flag);
      
      for(j=0;j<height;j++)
	{
          redata[j*width + i]= re[j];
	  imdata[j*width + i]= im[j];
	}
    }

  free(im);  free(re);

}


int calc_power_of_two(int number)
{
  int power;

  power =0;
  while (number !=1)
    {
      power++; number=number>>1;
    }
  return power;
}

void make_initial_data(double *re_part, double *im_part,
		       int num_of_data, int power)
{
  int i,j,ptr,offset,new_ptr,dft;
  double *rebuf, *imbuf;

  dft = num_of_data;
  rebuf = malloc(num_of_data * sizeof(double));
  imbuf = malloc(num_of_data * sizeof(double));

  for(i=1;i<power;i++){
    new_ptr = 0; offset = 0;
    while(new_ptr < num_of_data){
      ptr = 0;
      while (ptr<dft){
	rebuf[new_ptr]= *(re_part+offset+ptr);
	imbuf[new_ptr]= *(im_part+offset+ptr);
	new_ptr++;ptr+=2;
	if(ptr==dft)ptr =1;
      }
      offset +=dft;
    }
    for(j=0;j<num_of_data;j++){
      *(re_part+j)=rebuf[j];
      *(im_part+j)=imbuf[j];
    }
    dft /= 2;
  }
  free(rebuf);  free(imbuf);
}

void FFT1(double *re_part, double *im_part,
	  int num_of_data, int flag)
{
  int i,j,k,power,dft,half,num_of_dft,num_out,num_in1,num_in2;
  double unit_angle,step_angle,rebuf,imbuf,angle;

  static double *re_part_new,*im_part_new;

  re_part_new = malloc(num_of_data * sizeof(double));
  im_part_new = malloc(num_of_data * sizeof(double));

  if(flag == -1){
    for(i=0;i<num_of_data;i++){
      re_part [i] = re_part[i]/num_of_data;
      im_part [i] = -*(im_part+i)/num_of_data;
    }
  }

  power = calc_power_of_two(num_of_data);

  make_initial_data(re_part,im_part,num_of_data,power);
  unit_angle = 2.0 * PI / num_of_data;
  dft = 2;

  for(i=0;i<power;i++){
    num_of_dft = num_of_data/ dft;
    step_angle = unit_angle * num_of_dft;
    half = dft /2 ;

    for(j=0;j<num_of_dft;j++){
      angle = 0.0;

      for(k=0;k<dft;k++){
	num_out = j * dft +k;
	if(k<half)
	  {
	    num_in1 = num_out;
	    num_in2 = num_in1+half;
	  }
	else
	  {
	    num_in2 = num_out;
	    num_in1 = num_out - half;
	  }

	rebuf = *(re_part + num_in2);
	imbuf = *(im_part + num_in2);

	re_part_new[num_out] = 
	  *(re_part + num_in1) + rebuf * cos(angle) + imbuf * sin(angle);

	im_part_new[num_out] = 
	  *(im_part + num_in1) + imbuf * cos(angle) - rebuf * sin(angle);

	angle = angle + step_angle;
      }
    }

    for(j=0;j<num_of_data;j++)
      {
	*(re_part + j) = re_part_new[j];
	*(im_part + j) = im_part_new[j];
      }
    dft = dft * 2;
  }

  if(flag==-1)
    {
      for(j=0;j<num_of_data;j++){
	*(im_part+j)= - *(im_part + j);
      }
    }

  /********$B$3$N4X?t$O2?EY$b8F$P$l$k$N$G(B,
    free$B$7$F$*$/$3$H(B
    */
  free(im_part_new);
  free(re_part_new);
}