FFT(Fast Fourier Trasform)フーリエ変換の計算pie12.c

1200の自乗は1440000の計算


akira@dynabookAZ:~/pie$ gcc -o pie12 pie12.c -lm
akira@dynabookAZ:~/pie$ ./pie12
キーを押してください
0 0.000000 0.000000
1 0.000000 0.000000
2 2.000000 0.000000
3 1.000000 0.000000
4 0.000000 0.000000
5 0.000000 0.000000
6 0.000000 0.000000
7 0.000000 0.000000
8 0.000000 0.000000
9 0.000000 0.000000
10 0.000000 0.000000
11 0.000000 0.000000
12 0.000000 0.000000
13 0.000000 0.000000
14 0.000000 0.000000
15 0.000000 0.000000

0 0.000000 0.000000
1 0.000000 0.000000
2 0.000000 0.000000
3 -0.000000 0.000000
4 4.000000 0.000000
5 4.000000 -0.000000
6 1.000000 0.000000
7 0.000000 0.000000
8 0.000000 0.000000
9 -0.000000 -0.000000
10 -0.000000 0.000000
11 -0.000000 0.000000
12 0.000000 0.000000
13 -0.000000 0.000000
14 -0.000000 0.000000
15 0.000000 0.000000
time=0.047471
akira@dynabookAZ:~/pie$

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <termios.h>
#include <term.h>
#include <curses.h>
#include <math.h>
#include <time.h>

#define  P   4
#define  N   16

double f1[N+4],f2[N+4];

void  fft(double ff1[N+4],double ff2[N+4],int vec);
void a2(  double ff1[N+4],double ff2[N+4]);

static struct termios initial_settings, new_settings;
static int peek_character = -1;

void init_keyboard();
void close_keyboard();
int kbhit();
int readch();
double gettime();

double st1,pt1;
float tt1=0.0;
double w[N+N+4];
int n2=N>>1,n4=N>>2,n8=N>>3;

void  fft(double ff1[N+4],double ff2[N+4],int vec){
   int i,j,k,ik,h,d,k2;
   double t,s,c,dx,dy;

   j=0;
   for(i=0;i<=N-2;i++){
      if(i<j){
         t=ff1[j]; ff1[j]=ff1[i]; ff1[i]=t;
         t=ff2[j]; ff2[j]=ff2[i]; ff2[i]=t;
      }
      k=n2;
      while(k<=j){
         j-=k;k>>=1;
      }
      j+=k;
   }


   for(k=1;k<N; k=k2){
      h=0; k2=k+k; d=N/k2;
      for(j=0;j<k;j++){
         c=w[h+n4];
         if(vec!=-1) s=w[h];
         else              s=-w[h];
         for(i=j;i<N;i+=k2){
            ik=i+k;
            dx=s*ff2[ik]+c*ff1[ik];
            dy=c*ff2[ik]-s*ff1[ik];
            ff1[ik]=ff1[i]-dx; ff1[i]+=dx;
            ff2[ik]=ff2[i]-dy; ff2[i]+=dy;
         }
         h+=d;
      }
   }
   if(vec==-1){
      for(j=0;j<N;j++){
         ff1[j]/=N;
         ff2[j]/=N;
      }
   }
}

void  a2(double ff1[N+4],double ff2[N+4]){
   double gg1,gg2;
   int j;

   for(j=0;j<N;j++){
      gg1=ff1[j]*ff1[j]-ff2[j]*ff2[j];
      gg2=2*ff1[j]*ff2[j];
      ff1[j]=gg1;
      ff2[j]=gg2;
   }
}

int main()
{
   int ch=0,i,j;
   long cnt;
   double c,s,dc,ds,t;

   init_keyboard();
   printf("キーを押してください\n");

   st1=gettime();

      f1[2]=2.0;
      f1[3]=1.0;
      for(j=0;j<N;j++) printf("%d %6.6f %6.6f\n",j,f1[j],f2[j]);
      printf("\n");

   t=sin(M_PI/N);
   dc=2*t*t; ds=sqrt(dc*(2-dc));
   t=2*dc; c=w[n4]=1; s=w[0]=0;
   for(i=1;i<n8;i++){
      c-=dc; dc+=t*c;
      s+=ds; ds-=t*s;
      w[i]=s;
      w[n4-i]=c;
   }
   if(n8!=0) w[n8]=sqrt(0.5);
   for(i=0;i<n4;i++)    w[n2-i]=w[i];
   for(i=0;i<n2+n4;i++) w[i+n2]=-w[i];

   cnt=0;
   while(ch != 'q'&&cnt<100) {
      for(j=0;j<N;j++){
         f1[j]=f2[j]=0.0;
      }
      f1[2]=2.0;
      f1[3]=1.0;

      fft(f1,f2,1);
      a2( f1,f2);
      fft(f1,f2,-1);

      if(kbhit()) {
         ch = readch();
         printf("you hit %X : %c\n",ch,ch);
      }
      cnt++;
   }
   for(j=0;j<N;j++) printf("%d %6.6f %6.6f\n",j,f1[j],f2[j]);

   pt1=(gettime()-st1)+tt1;
   printf("time=%6.6f\n",pt1);

    close_keyboard();
}

void init_keyboard()
{
    tcgetattr(0, &initial_settings);
    new_settings = initial_settings;
    new_settings.c_lflag &= ~ICANON;
    new_settings.c_lflag &= ~ECHO;
    new_settings.c_lflag &= ~ISIG;
    new_settings.c_cc[VMIN] = 0;
    new_settings.c_cc[VTIME] = 0;
    tcsetattr(0, TCSANOW, &initial_settings);
}

void close_keyboard()
{
    tcsetattr(0, TCSANOW, &initial_settings);
}

int kbhit()
{
    char ch;
    int nread;

    if(peek_character != -1)
        return 1;
    new_settings.c_cc[VMIN]=0;
    tcsetattr(0, TCSANOW, &new_settings);
    nread = read(0, &ch, 1);
    new_settings.c_cc[VMIN]=1;
    tcsetattr(0, TCSANOW, &new_settings);

    if(nread == 1) {
        peek_character = ch;
        return 1;
    }
    return 0;
}


int readch()
{
    char ch;

    if(peek_character != -1) {
        ch = peek_character;
        peek_character = -1;
        return ch;
    }
    read(0, &ch, 1);
    return ch;
}


double gettime(){
   struct timeval tv;
   gettimeofday(&tv,NULL);
   return tv.tv_sec+(double)tv.tv_usec*1e-6;
}