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
113 TCP/IP State Transition - TCP 스택 포팅 시 참조 file digipine 2017.11.02 197130
112 언어 IDE 별로 git ignore 파일을 자동으로 만들어 주는 사이트 엉뚱도마뱀 2018.12.17 126067
111 What is Android Repo? and Manual, Download file digipine 2017.11.02 100764
110 [Linux] ubuntu 16.04에 QT Creator 설치하기 digipine 2017.11.02 24872
109 [Swift, MacOS] 맥 한글 파일명이 윈도우에서 자소 분리되는 현상 해결, NFD, NFC 엉뚱도마뱀 2018.12.11 20544
108 Linux init.d 에서 등록하기. 부팅 시 자동실행 설정 digipine 2017.11.03 13474
» FFT (Fast Fourier Transform) 고속 푸리에 변환 예제 소스 digipine 2017.10.29 12841
106 Photoshop CC 2018 한글 영문 변환 언어팩, 포토샵 언어변경 file 엉뚱도마뱀 2018.07.04 8866
105 Phabricator 설치 가이드 우분투 16.04 기준 digipine 2017.11.02 6721
104 ATmega8 MCU 간의 TWI 기능을 이용한 I2C 통신 digipine 2017.11.02 6511
103 WinPCap과 Ethereal, Wireshark 을 이용한 스니핑(Sniffing) digipine 2017.10.29 5819
102 난수발생기 개론 엉뚱도마뱀 2017.11.22 4787
101 비밀번호 해쉬에 Salt(소금) 첨가하기 file 엉뚱도마뱀 2017.11.23 4669
100 공짜 무료 C/C++ 컴파일러들 file digipine 2017.10.28 4666
99 [iOS] Bluetooth로 App을 백그라운드 모드로 실행는 방법 lizard2019 2020.02.11 3949
98 [Linux, OSX] pfctl - Packet FIlter Control 사용법 digipine 2017.11.02 3767
97 OpenAL 사용법 정리 1 digipine 2017.11.01 3733
96 WIN CE, GPS - NMEA protocol - GPS Virtual Driver digipine 2017.10.28 3310
95 Xcode 없이 맥에 '명령어 라인 도구(Command Line Tools)'를 설치하는 방법 엉뚱도마뱀 2018.12.26 3279
94 LibVLC 미디어 재생기 프로그래밍 방법 C++, QT 엉뚱도마뱀 2018.04.20 2927
Board Pagination Prev 1 2 3 4 5 6 Next
/ 6