logo

English

이곳의 프로그래밍관련 정보와 소스는 마음대로 활용하셔도 좋습니다. 다만 쓰시기 전에 통보 정도는 해주시는 것이 예의 일것 같습니다. 질문이나 오류 수정은 siseong@gmail.com 으로 주세요. 감사합니다.

FFT (Fast Fourier Transform) 고속 푸리에 변환 예제 소스

by digipine posted Oct 29, 2017
?

Shortcut

PrevPrev Article

NextNext Article

Larger Font Smaller Font Up Down Go comment Print
?

Shortcut

PrevPrev Article

NextNext Article

Larger Font Smaller Font Up Down Go comment Print

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

#define MAX(a,b)    (((a) > (b)) ? (a) : (b))
#define PI  3.141592

typedef struct {
    double real, imag;
} COMPLEX;

   FILE *fp;
 
/* function prototypes for DSP disk file functions */
   FILE *open_read(char *);
   FILE *open_write(char *);
   void read_record(FILE *,double *,int);
   void write_record(double *,int);
   double *read_double_record(FILE *);

/* function prototypes for dft and inverse dft functions */
   void fft(COMPLEX *,int);
   int log2(unsigned int);

/************************************************************************
IFFTTEST.C - Demonstrate and test FFT and Inverse FFT functions
Requires time domain data file in DSP_FILE format.
Generates DSP_FILE format files for spectral magnitude
and reconstructed time domain data.
************************************************************************/

void main(void)
{
    int          i, length, fft_length, m;
    double       tempflt;
    double       *signal,*log_mag;
    COMPLEX      *samp;

// Read the input data file from the dsp format.  데이타파일 읽는 명령
    fp = open_read("test.txt");
    length = 128;
    signal = read_double_record(fp);

// determine fft size and allocate the complex array - fft size와 배열
    m = log2(length);
    fft_length = 1 << m;   //bit shift
    samp = (COMPLEX *) calloc(fft_length, sizeof(COMPLEX));

// copy input signal to complex array and do the fft
    for (i=0; i<length; i++)
  samp[i].real = signal[i];
    fft(samp, m);

/* find log magnitude and store for output */
    log_mag = (double *)calloc(fft_length, sizeof(double));
    for (i=0; i< fft_length; i++)
 {
        tempflt  = samp[i].real * samp[i].real;
        tempflt += samp[i].imag * samp[i].imag;
        tempflt =  sqrt(tempflt);
        log_mag[i] = 10 * log10(MAX(tempflt, 1.e-14));
    }
    fp = open_write("test_out.txt");  //결과 저장 텍스트 //
    write_record(log_mag, length);
}

/**************************************************************************
DISKIO.C - Source code for DSP data format read and write functions
    open_read          open DSP data file to be read
    open_write         create header and open DSP data file for write
    read_record        read one record
    write_record       write one record
    read_float_record  read one record and convert to float array
*************************************************************************/

/***************************************************************************************
open_read - open a DSP data file for read
Returns a pointer to a DSP_FILE structure allocated by the function and opens file_name.
Allocation errors or improper type causes a call to exit(1).
A bad file_name returns a NULL pointer.
DSP_FILE *open_read(char *file_name)
***************************************************************************************/

FILE *open_read(char *file_name) /* file name string */
{
/* allocate the DSP data file structure */

    fp = (FILE *) malloc(sizeof(FILE));

/* open file for text read and update */
    fp = fopen(file_name,"r");
    return(fp);
}

/**************************************************************************
open_write - open a DSP data file for write
Returns a pointer to a DSP_FILE structure allocated by the function.
Allocation errors or improper type causes a call to exit(1).
A bad file name returns a NULL pointer.
DSP_FILE *open_write(char *file_name, int records,int rec_len)
    file_name       pointer to file name string
    records         number of records of data to be written
    rec_len         number of elements in each record
*************************************************************************/

FILE *open_write(char * file_name)
{
/* allocate the DSP data file structure */
    fp = (FILE *) malloc(sizeof(FILE));

/* open file for text write and update*/
    fp = fopen(file_name,"wt");
    return(fp);
}

/********************************************************************
read_record - read one record of DSP data file
Exits if a read error occurs or if the DSP_FILE structure is invalid.   
Void read_record(FILE *fp , double *input , int length)
********************************************************************/

void read_record(FILE *fp ,double * input , int length)
{
    int status , i;
 for(i=0 ; i < length ; i++)
 {
  status = fscanf(fp,"%lf\n", &input[length + i]);
    }
}

/************************************************************************************************************************
read_double_record - read one record of DSP data file and convert to float array of values.
Returns a pointer to the beginning of the allocated float array of values representing the record read from the DSP_FILE.
Exits if a read or allocation error occurs.
float *read_double_record(DSP_FILE *dsp_info)
************************************************************************************************************************/

double *read_double_record(FILE * fp)
{
    static double *buf;            /* input buffer to read data in */
    double *out;                    /* return output pointer */
    double *out_ptr;
 double *d_ptr;
    int i,length,length_in;

  length = 128;
    length_in = 256;
 
    buf = (double *) calloc(length_in,sizeof(double));
    out = (double *) calloc(length,sizeof(double));

/* read the record into buf */
    read_record(fp , buf , length);

/* perform conversion to floating point */
     out_ptr = out; 
     d_ptr = buf + 128;
     for(i = 0 ; i < length ; i++)
         *out_ptr++ = (double)(*d_ptr++);
 
  return(out);                 /* return converted pointer */
}

/**************************************************************************
write_record - write one record of DSP_FILE data
Exits if write error occurs or if the DSP_FILE structure is invalid.
void write_record(char *ptr,DSP_FILE *dsp_info, int length)
    ptr        pointer to data to write to disk (type in dsp_info)
    dsp_info   pointer to DSP data file structure   
*************************************************************************/

void write_record(double *ptr,int length)
{
    int status, i;
 for(i=0 ; i < length ; i++)
 {
  status = fprintf(fp,"%lf\n", *ptr++);
    }
}

/**********************************************************************
DFT.C - SOURCE CODE FOR DISCRETE FOURIER TRANSFORM FUNCTIONS
fft     In-place radix 2 decimation in time FFT
log2    Base 2 logarithm
***********************************************************************/

/**********************************************************************
fft - In-place radix 2 decimation in frequency FFT
Requires pointer to complex array, x and power of 2 size of FFT, m
(size of FFT = 2^m).  Places FFT output on top of input COMPLEX array.
void fft(COMPLEX *x, int m)
***********************************************************************/

void fft(COMPLEX *x, int m)

  static COMPLEX *w;           // used to store the w complex array
    static int mstore = 0;       // stores m for future reference
    static int n = 1;            // length of fft stored for future

    COMPLEX u, temp, tm;
    COMPLEX *xi, *xip, *xj, *wptr;

    int i, j, k, l, le, windex;

    double arg, w_real, w_imag, wrecur_real, wrecur_imag, wtemp_real;

    if(m != mstore)
 {
  // free previously allocated storage and set new m
        if(mstore != 0) free(w);
        mstore = m;
        if(m == 0) return;       // if m=0 then done

  // n = 2^m = fft length
        n = 1 << m;  
        le = n/2;  //difference between the upper and lower leg indices

  // allocate the storage for w
        w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));

  // calculate the w values recursively 계산
        arg = PI/le;         //  PI/le calculation
        wrecur_real = w_real = cos(arg);
        wrecur_imag = w_imag = -sin(arg);
        xj = w;
        for (j = 1 ; j < le ; j++)
  {
            xj->real = (double)wrecur_real;
            xj->imag = (double)wrecur_imag;
            xj++;
            wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
            wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
            wrecur_real = wtemp_real;
        }
    }
 // start fft 푸리에 변환 시작
    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++)
 {
        le = le/2;
  // first iteration with no multiplies
        for(i = 0 ; i < n ; i = i + 2*le)
  {
            xi = x + i;
            xip = xi + le;
            temp.real = xi->real + xip->real;
            temp.imag = xi->imag + xip->imag;
            xip->real = xi->real - xip->real;
            xip->imag = xi->imag - xip->imag;
            *xi = temp;
        }
  // remaining iterations use stored w
        wptr = w + windex - 1;
        for (j = 1 ; j < le ; j++)
  {
            u = *wptr;
            for (i = j ; i < n ; i = i + 2*le)
   {
                xi = x + i;
                xip = xi + le;
                temp.real = xi->real + xip->real;
                temp.imag = xi->imag + xip->imag;
                tm.real = xi->real - xip->real;
                tm.imag = xi->imag - xip->imag;            
                xip->real = tm.real*u.real - tm.imag*u.imag;
                xip->imag = tm.real*u.imag + tm.imag*u.real;
                *xi = temp;
            }
            wptr = wptr + windex;
        }
        windex = 2*windex;
    }
 // rearrange data by bit reversing
    j = 0;
    for (i = 1 ; i < (n-1) ; i++)
 {
        k = n/2;
        while(k <= j)
  {
            j = j - k;
            k = k/2;
        }
        j = j + k;
        if (i < j)
  {
            xi = x + i;
            xj = x + j;
            temp = *xj;
            *xj = *xi;
            *xi = temp;
        }
    }
}

/**************************************************************************
log2 - base 2 logarithm
Returns base 2 log such that i = 2^ans where ans = log2(i).
if log2(i) is between two values, the larger is returned.
int log2(unsigned int x)
*************************************************************************/
int log2(unsigned int x)
{
    unsigned int mask,i;

    if(x == 0) return(-1);      // zero is an error, return -1
    x--;                        // get the max index, x-1
    for(mask = 1 , i = 0 ; ; mask *= 2 , i++)
 {
        if(x == 0) return(i);   // return log2 if all zero
        x = x & (~mask);        // AND off a bit
    }

}


http://starpotato.springnote.com/pages/3236324 권영준님의 노트에서...

Jan.06, 2009 : FFT의 개괄적인 이해는 2d fast fourier transform 을 참조하세요

FFT 소스를 요청하는 분들이 많이 계셔서 소스코드를 공개합니다. 이게 기본틀은 어디서 코드를 분석한 것이라서 거의 똑같다고 보시면 됩니다. 간단한 주석을 달았구요. 도움이 되시길 빕니다.

          /* File            :       struct.h
           * Creator      :       BlackEngine
           * Date         :       2006/03/28
           * Version      :       0.0.1 ( created )
           * Descript     :       Definition of Complex Number , and etc
           */
           
          typedef struct _COMPLEX {
                  double     real;
                  double     imag;
          } COMPLEX, *pCOMPLEX;
          
          /* File            :       fft.h
           * Creator      :       BlackEngine
           * Date         :       2006/03/28
           * Version      :       0.0.1 ( created )
           * Descript     :       Perform the FFT Algorithm
           */
          #ifndef _BLACKENGINE_FFT_H_
          #define _BLACKENGINE_FFT_H_
           
          #include "struct.h"
          /* TODO : at 03/29 , decision the return values of FFT Function. */
          int fft_2d(COMPLEX ** ,int ,int ,int );
          int fft(int ,int ,double *,double *);
          #endif /* ifndef _BLACKENGINE_FFT_H_ */
          
          /* File            :       fft.c
           * Creator      :       BlackEngine
           * Date         :       2006/03/28
           * Version      :       0.0.1 ( created )
           * Descript     :       Perform the FFT Algorithm ( 2-D Image FFT )
           */
          #include 
          #include 
          #include
           
          #include "fft.h"
           
          #define TRUE 1
          #define FALSE 0
          /*#define NULL 0*/
           
          /* Powerof2 function : compare the integer with 2^m
           * and find maximum m , that is 2^m &lt;= nx
           * twopm = 2^m
           */
          int Powerof2(int nx, int* m, int* twopm) {
                  int pwr;
                  *m = 0;
                  for ( pwr = 1 ; pwr &lt; nx ; pwr = pwr*2) {
                          *m = *m + 1;
                  }
                  *twopm = pwr;
                  return(TRUE);
          }
           
          int fft_2d(COMPLEX **c, int nx, int ny, int dir) {
                  int i,j;
                  int m,twopm;
                  double *real,*imag;
           
                  /* Transform the rows */
                  real = (double *)malloc(nx * sizeof(double));
                  imag = (double *)malloc(nx * sizeof(double));
                  if (real == NULL || imag == NULL)
                          return(FALSE);
                  if (!Powerof2(nx,&amp;m,&amp;twopm) || twopm != nx)
                          return(FALSE);
                  for (j=0;j&lt;&lt;= m;
           
                  /* N의 절반값을 half에 저장한다. */
                  half = N &gt;&gt; 1;
           
                  /* Bit Reverse를 하는 과정 */
                  for ( i = 0 , j = 0 ; i &lt; N - 1 ; i++ ) {
                          if ( i &lt; j )
                          {
                                  swap_x = x[i];
                                  swap_y = y[i];
                                  x[i] = x[j];
                                  y[i] = y[j];
                                  x[j] = swap_x;
                                  y[j] = swap_y;
                          }
                          step = half;
           
                          while ( step &lt;= j )
                          {
                                  j = j - k;
                                  k &gt;&gt;= 1;
                          }
                          j = j + step;
                  }
           
                  /* 초기의 값을 입력하는 부분,
                   * factor의 역할은 각 단계(Stage)에서 omega의 값을 변경하기 위한
                   * 기본 값이다. 즉 omega N= 4일 경우에는 각도가 pi/2 만큼씩 변하므로
                   * exp(-j pi/2) 의 값을 가지고 있는 것이 factor이다
                   * factor는 각 stage가 올라갈 때 마다 변경되어야 한다.
                   * 그 각도가 절반으로 줄어야 하므로
                   * cos(t) + j sin(t) 가 cos(t/2) + jsin(t/2) 로 되어야 한다.
                   * cos(t) = 2cos^2(t/2) -1 이므로 cos(t/2) = sqrt((1+cos(t))/2)이고
                   * sin(t/2) = sqrt(1-cos^2(t/2))이므로 sqrt((1-cos(t))/2) 이다.
                   */
                  factor_x = -1.0;
                  factor_y = 0;
                  butterfly = 1;
           
                  for ( stage = 0 ; stage &lt; m ; stage++ )
                  {
                          /* step 은 butterfly가 두개의 샘플을 가지고 작업하므로
                           * 현재 위치에서 몇번째 샘플을 가지고 해야 하는 가를
                           * 가리키는 변수이다.
                           * butterfly는 같은 omega 값에의해 반복되는 butterfly가
                           * 존재하므로 그 위치를 가리키기 위한 값이다.
                           */
                          step = butterfly;
                          butterfly &lt;&lt;= 1;
                          /* 각 단계에서 omega는 항상 1+j0으로 시작한다. */
                          omega_x = 1.0;
                          omega_y = 0;
                          for ( numofbutter = 0 ; numofbutter &lt; step ; numofbutter++)
                          {
                                  for( i = numofbutter; i &lt; N; i = i + butterfly)
                                  {
                                          /* omega_N(k) 일때 omega_N(k+N/2) = -omega_N(k)
                                           * 의 성질을 이용한다. 그래서 두번째 샘플에 미리
                                           * omega_N(k)를 곱한 후 결과의 처음 샘플에는
                                           * 그냥 더하고 두번째 샘플에는 빼준다(-이므로)
                                           */
                                          swap_x = x[i+step]* omega_x - y[i+step]*omega_y;
                                          swap_y = y[i+step]*omega_x + x[i+step]*omega_y;
           
                                          x[i+step] = x[i] - swap_x;
                                          y[i+step] = y[i] - swap_y;
           
                                          x[i] = x[i] + swap_x;
                                          y[i] = y[i] + swap_y;
           
                                  }
                                  /* 공통된 omega를 사용하는 샘플을 다 곱한 후
                                   * omega를 factor의 각도만큼 증가시켜야 한다.
                                   * omega(k+1) 의 값은 omega와 factor의 곱을 통해
                                   * 구할 수 있다.
                                   */
                                  swap_x = omega_x * factor_x - omega_y * factor_y;
                                  swap_y = omega_y * factor_x + omega_x * factor_y;
                                  omega_x = swap_x;
                                  omega_y = swap_y;
           
                          }
                          /* 한 단계가 끝났으므로 factor를 기존 값에서 각도가
                           * 절반인 값으로 구해야 한다.
                           */
                          factor_y = sqrt((1.0+factor_x)/2.0);
                          if( dir == 1)
                                  factor_y = -factor_y;
                          factor_x = sqrt((1.0-factor_x)/2.0);
                  }
                  /* Inverse FFT 의 경우에는 sampling을 해줘야 한다.
                   * 다른 소스에서는 FFT의 값이 크게 나오므로 fft시에 sampling을
                   * 하고 inverse에서 안하는 경우도 있다 ( 상관없음 )
                   */
                  if( dir == -1 )
                  {
                          for ( i = 0 ; i &lt; N ; i++)
                          {
                                  x[i] /= N;
                                  y[i] /= N;
                          }
                  }
           
                  return TRUE;
          }
          


=======dft.h============

typedef struct {
	double real, imag;
} COMPLEX;
#define PI 3.141592

void fft(COMPLEX *,int);
void ifft(COMPLEX *,int);
void dft(COMPLEX *,COMPLEX *,int);
void idft(COMPLEX *,COMPLEX *,int);
void rfft(double *,COMPLEX *,int);
void ham(COMPLEX *,int);
void han(COMPLEX *,int);
void triang(COMPLEX *,int);
void black(COMPLEX *,int);
void harris(COMPLEX *,int);
int log2(unsigned int);


=========dft.c==========

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

void fft(COMPLEX *x, int m)
{ 
	static COMPLEX *w; 
	static int mstore = 0;
	static int n = 1;
	COMPLEX u,temp,tm;
	COMPLEX *xi,*xip,*xj,*wptr;
	int i,j,k,l,le,windex;
	double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
	if(m != mstore) {

		if(mstore != 0) free(w);
		mstore = m;
		if(m == 0) return; 
		n = 1 << m; 
		le = n/2; 

		w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
		if(!w) {
			printf("\nUnable to allocate complex W array\n");
			exit(1);
		}

		arg = PI/le; 
		wrecur_real = w_real = cos(arg);
		wrecur_imag = w_imag = -sin(arg);
		xj = w;
		for (j = 1 ; j < le ; j++) {
			xj->real = (double)wrecur_real;
			xj->imag = (double)wrecur_imag;
			xj++;
			wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
			wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
			wrecur_real = wtemp_real;
		}
	}

	le = n;
	windex = 1;
	for (l = 0 ; l < m ; l++) {
		le = le/2;

		for(i = 0 ; i < n ; i = i + 2*le) {
			xi = x + i;
			xip = xi + le;
			temp.real = xi->real + xip->real;
			temp.imag = xi->imag + xip->imag;
			xip->real = xi->real - xip->real;
			xip->imag = xi->imag - xip->imag;
			*xi = temp;
		}


		wptr = w + windex - 1;
		for (j = 1 ; j < le ; j++) {
			u = *wptr;
			for (i = j ; i < n ; i = i + 2*le) {
				xi = x + i;
				xip = xi + le;
				temp.real = xi->real + xip->real;
				temp.imag = xi->imag + xip->imag;
				tm.real = xi->real - xip->real;
				tm.imag = xi->imag - xip->imag; 
				xip->real = tm.real*u.real - tm.imag*u.imag;
				xip->imag = tm.real*u.imag + tm.imag*u.real;
				*xi = temp;
			}
			wptr = wptr + windex;
		}
		windex = 2*windex;
	} 

	j = 0;
	for (i = 1 ; i < (n-1) ; i++) {
		k = n/2;
		while(k <= j) {
			j = j - k;
			k = k/2;
		}
		j = j + k;
		if (i < j) {
			xi = x + i;
			xj = x + j;
			temp = *xj;
			*xj = *xi;
			*xi = temp;
		}
	}
}

main(){
	int i;
	COMPLEX com[1];
	fft(com,1);
	for(i=0;i<sizeof(com);i++)
		printf("%f\n",com[i].real);
}
 
TAG •

List of Articles
No. Subject Author Date Views
112 공짜 무료 C/C++ 컴파일러들 file digipine 2017.10.28 4594
111 UML - 기본편 ( 기본 표기 형식 및 관계 표현법 ) digipine 2017.10.28 1056
110 프로세스 능력 성숙도 모델(CMMI)의 적용 digipine 2017.10.28 893
109 ALM의 등장 배경, 오해와 진실 digipine 2017.10.28 522
108 CMM / CMMI 란 무엇인가? digipine 2017.10.28 2739
107 Thread간 동기화관리자 digipine 2017.10.28 759
106 Windows API - 안전한 문자열 함수들 digipine 2017.10.28 537
105 WIN CE C++ 시리얼 제어 방법 digipine 2017.10.28 908
104 WIN CE, GPS - NMEA protocol - GPS Virtual Driver digipine 2017.10.28 3171
103 NMEA-0183을 이용한 GPS 애플리케이션 제작 digipine 2017.10.28 728
102 Windows에서 SVN 용 폴더 한 번에 삭제하기 digipine 2017.10.29 771
101 WinPCap과 Ethereal, Wireshark 을 이용한 스니핑(Sniffing) digipine 2017.10.29 5685
» FFT (Fast Fourier Transform) 고속 푸리에 변환 예제 소스 digipine 2017.10.29 12738
99 Linux(리눅스) 파일 시스템 정리 (ext, ext2, ext3, ext4) digipine 2017.10.29 1273
98 Ubuntu Server OS 한국어 모드로 설치 후 서버 콘솔에서 한글 깨짐 해결 방법 digipine 2017.10.31 2454
97 QCAD 1.1 한국 최초의 2D CAD 소프트웨어 file digipine 2017.11.01 1051
96 OpenAL 사용법 정리 1 digipine 2017.11.01 3603
95 WPA_SUPPLICANT 빌드 방법 digipine 2017.11.01 540
94 WPA_SUPPLICANT 빌드 방법 digipine 2017.11.01 723
93 Iconv 사용법 소스 digipine 2017.11.01 1536
Board Pagination Prev 1 2 3 4 5 6 Next
/ 6