//////// CZT and ICZT: Chirp Z-transform Full Source Code For Image Data /////////
//////// Smart Fast Fourier Transform for any size image, Full Source Code /////////
// CZT_TwoDimensional.c (+ stb_image.h), ONLY single file source code. //
// Everything (DFT, FFT, CZT, etc.) is implemented within it, //
// So NO external libraries are used except for stb_image.h //
// for reading image file. //
// //
// Multi-threaded version of Fourier Transform source code for any size of 2D data. //
// //
// gcc -fopenmp -O3 -Wall -W -o CZT_TwoDimensional CZT_TwoDimensional.c -lm //
// //
// //
// You have the right to: You are free to use this source code and modify it //
// for your own use without any restrictions. In addition, we do not provide //
// any warranty or support under any circumstances. //
// //
// //
//////// https://github.com/nothings/stb/blob/master/stb_image.h
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include <complex.h>
#undef PI
#define PI 3.1415926535897932384626433832795L // 円周率(π)
#define Rotor(a) (double complex)( cos(a) + I*sin(a) )
#define COMPLEX(a,b) (double complex)( (a) + I*(b) )
enum { DFT, IDFT };
enum { FFT=DFT, IFFT=IDFT };
enum { CZT=DFT, ICZT=IDFT };
#pragma pack(push,1)
typedef struct PIXEL { unsigned char B, G, R, A; } PIXEL ;
#pragma pack(pop)
typedef struct COMPLEX { double r, i; } COMPLEX;
#pragma pack(push,2)
typedef struct HEADER { // BMP => bfType==0x4D42 (0x42,0x4D) little endian
unsigned short bfType; // identification 0x42='B', 0x4D='M'
unsigned int bfSize; // ??
unsigned short bfReserved1; // 0
unsigned short bfReserved2; // 0
unsigned int bfOffBits; // 54 sizeof(HEADER)
////////////////////////////////////////////////////////////////////////////
unsigned int biSize; // 40 sizeof(biSize ～ biClrImportant)
signed int biWidth; // lx
signed int biHeight; // ly
unsigned short biPlanes; // 1
unsigned short biBitCount; // 32
unsigned int biCompression; // 0
unsigned int biSizeImage; // ??
signed int biXPelsPerMeter; // 0
signed int biYPelsPerMeter; // 0
unsigned int biClrUsed; // 0
unsigned int biClrImportant; // 0
////////////////////////////////////////////////////////////////////////////
unsigned char lut[];
} HEADER ;
#pragma pack(pop)
#define BMPID 0x4D42
#define BMPIFHSIZE 40
#define saveArray2bmp32(name,a,w,h) \
saveDataAsBMP32((name),(&(a)[0][0]),(w),(h))
#define IF if
#define EI else if
#define EE else
int dft=0, verbose=0;
void **alloc2D( void *data, int s, int lx, int ly ){
long long i, w=lx, h=ly; if( !data || s<1 || lx<1 || ly<1 ) return NULL;
void **d, *base=data ;
if( (d=(void **)malloc(h*sizeof(void *))) == NULL ) return NULL ;
for( i=0; i<ly; i++ ) d[i] = (void *)((unsigned char *)base + w*s*i) ;
return d ;
}
void **alloc2Darray( int s, int lx, int ly ){
long long i, w=lx, h=ly; if(s<1 || lx<1 || ly<1) return NULL ;
void **d, *base ;
if( (d=(void **)calloc(1,h*sizeof(void *)+(w*h*s))) == NULL ) return NULL ;
base = (unsigned char *)d + h*sizeof(void *) ;
for( i=0; i<ly; i++ ) d[i] = (void *)((unsigned char *)base + w*s*i) ;
return d ;
}
int saveDataAsBMP32( char *fname, void *b, int lx, int ly ){
HEADER h = {BMPID, lx*abs(ly)*sizeof(PIXEL)+sizeof(HEADER), 0, 0, sizeof(HEADER),
BMPIFHSIZE, lx, ly, 1, 32, 0, lx*abs(ly)*sizeof(PIXEL), 0, 0, 0, 0};
FILE *fp=NULL ;
int st=-1;
if( (fp=fopen(fname, "wb")) == NULL ) {
fprintf( stderr, "File '%s': cannot create.\n", fname );
} else if ( fwrite( &h, sizeof(HEADER), 1, fp ) != 1 ) {
fprintf( stderr, "File '%s': Header write error.\n", fname );
} else if ( fwrite( b, h.biSizeImage, 1, fp ) != 1 ) {
fprintf( stderr, "File '%s': Data write error.\n", fname );
} else st=0;
if(fp)fclose(fp) ;
return st;
}
////////////////////////////////////////////////////////////////////////////////
////
//// 二次元の離散フーリエ変換
//// Two-dimensional discrete Fourier transform (DFT)
////
#include <omp.h> // gcc -fopenmp -O3 ( -O3: Strongly recommendation )
#define AppendCOMPLEX(a,b) { (a).r += (b).r, (a).i += (b).i ; }
#define ConjCOMPLEX(a) (COMPLEX){ (a).r, -((a).i) }
#define RotorCOMPLEX(a) (COMPLEX){ cos(a), sin(a) }
#define AddCOMPLEX(a,b) (COMPLEX){ (a).r + (b).r, (a).i + (b).i }
#define SubCOMPLEX(a,b) (COMPLEX){ (a).r - (b).r, (a).i - (b).i }
#define MulCOMPLEX(a,b) (COMPLEX){ (a).r*(b).r - (a).i*(b).i, (a).r*(b).i + (a).i*(b).r }
#define ImeMulCOMPLEX(a,b) (COMPLEX){ (a).r * (b), (a).i * (b) }
extern int verbose;
void dft_TwoDimensional( COMPLEX **d, int lx, int ly, int flag ) {
if(verbose)printf("execution: %s\t%d\n", __func__, flag );
int x, y, elements=lx;
double unit;
COMPLEX *rotor = (COMPLEX *)calloc( (lx>ly?lx:ly), sizeof(COMPLEX) );
if(!rotor) return;
elements = lx; unit = (flag==IDFT) ? 2*M_PI/elements : -2*M_PI/elements;
for( int e=0; e<elements; rotor[e]=RotorCOMPLEX(unit*e), e++ );
#pragma omp parallel for
for( y=0; y<ly; y++ ) { COMPLEX tmp[elements];//You need to be aware of the thread stack overflow.
for(int u=0; u<elements; u++ ){ int v; COMPLEX t = (COMPLEX){0,0};
for( v=0; v<elements; v++ ) AppendCOMPLEX( t, MulCOMPLEX(d[y][v],rotor[(u*v)%elements]) );
tmp[u] = t;
} double reciprocal = 1./elements; // tmp[u]/elements -> tmp[u]*reciprocal
if(flag==IDFT) for(int e=0; e<elements; e++ ) d[y][e] = ImeMulCOMPLEX(tmp[e],reciprocal);
else for(int e=0; e<elements; e++ ) d[y][e] = tmp[e];
} //////////////////////////////////////////////////////////////////////////
elements = ly; unit = (flag==IDFT) ? 2*M_PI/elements : -2*M_PI/elements;
if(lx!=ly){ for( int e=0; e<elements; rotor[e]=RotorCOMPLEX(unit*e), e++ ); }
#pragma omp parallel for
for( x=0; x<lx; x++ ) { COMPLEX tmp[elements];//You need to be aware of the thread stack overflow.
for(int u=0; u<elements; u++ ){ int v; COMPLEX t = (COMPLEX){0,0};
for( v=0; v<elements; v++ ) AppendCOMPLEX( t, MulCOMPLEX(d[v][x],rotor[(u*v)%elements]) );
tmp[u] = t;
} double reciprocal = 1./elements; // tmp[u]/elements -> tmp[u]*reciprocal
if(flag==IDFT) for(int e=0; e<elements; e++ ) d[e][x] = ImeMulCOMPLEX(tmp[e],reciprocal);
else for(int e=0; e<elements; e++ ) d[e][x] = tmp[e];
} //////////////////////////////////////////////////////////////////////////
free(rotor);
}
////////////////////////////////////////////////////////////////////////////////
// Twiddle factor table
typedef struct TWIDDLE_FACTOR { int elements, flag; double complex *table; }
TWIDDLE_FACTOR, TWIDDLE_TABLE ;
////////////////////////////////////////////////////////////////////////////////
////
//// CZT Chirp Z-transform, Bluestein's algorithm
////
//// Structure to store data for CZT Chirp Z-transform
//// samples : number of actual data
//// samples_out: Actual data count (for output, same as input one)
//// samples_fft: Number of data used in CZT conversion (power of 2)
//// flag_fft : flag indicating if "samples" is a power of 2
//// W[2] : For FFT/IFFT, Twiddle factor table
//// w : The number of weight data elements is (samples)
//// v : The number of impulse response data elements is (samples_fft)
////
typedef struct CZT_CONTROL {
int samples; // 標本点の数
int samples_out; // 出力する標本点の数
int samples_fft; // 2の整数乗値：flag_fft?samples:2*get_czt_size(samples)
int flag_fft; // 直接FFTを実行できるサイズかの判定
TWIDDLE_TABLE W[2];// 回転因子用テーブル
double complex *w; // 重みデータ 要素数は(samples)
double complex *v; // インパルス応答データ 要素数は(samples_fft)
} CZT_CONTROL;
////////////////////////////////////////////////////////////////////////////////
//// checks if the argument is a power of 2 ( Returns 1 if it is a power of 2 )
////
int ft_size_check( int n ){ if( n<4 ) return 0;
for(int e=n;e>1;e>>=1 ) if( e&1 ) return 0;
return 1;
}
////////////////////////////////////////////////////////////////////////////////
//// A function that returns the smallest power of 2 that is greater than or
//// equal to twice its argument
////
int get_czt_size( int n ){
if( ft_size_check(n) ) return n;
int a=1, b= n<<1;
for( a=1; a<b; a<<=1 );
return a;
}
////////////////////////////////////////////////////////////////////////////////
//// FFT routine for simple one-dimensional data
//// double complex *c: Pointer to complex array for input/output
//// elements : number of data in complex array
//// flag : Flag indicating forward transform or inverse transform
//// { FFT | IFFT }
////
void fft( double complex *c, int elements, int flag ){
if( !ft_size_check(elements) ) return; //2の冪(べき)のみ、power of two, ONLY
double complex *rotor = (double complex *)calloc( elements, sizeof(double complex) );
if(!rotor) return; // Processing is aborted because the twiddle factor table
// used for calculation could not be secured.
double unit = (flag==IFFT) ? 2*M_PI/elements : -2*M_PI/elements;;
for(int e=0; e<elements; rotor[e]=Rotor(unit*e), e++ );
for(int i=0, j=1; j<elements-1; j++ ) {// scrambler, for initial placement
for(int k=elements>>1; k>(i^=k); k>>=1 );
if( j<i ){double complex temp=c[j]; c[j]=c[i], c[i]=temp; } // exchange
}
for(int m, mh=1; (m=mh<<1)<=elements; mh=m ) {// Operations
for(int ir=0, i=0; i<elements; i+=m ) { // ir=0: Rotation angle is 0rad.
register double complex temp, W=rotor[ir]; // W: Rotation factor
for(int k=elements>>2; k>(ir^=k); k>>=1 );
for(int j=i; j<mh+i; j++ ) { int a = j, b = j + mh;
temp = c[a] - c[b]; // For butterfly operation
c[a] = c[a] + c[b]; // a ->- (a + b) ->- a
c[b] = temp * W; // b ->- (a-b)*W ->- b
} } } double reciprocal = 1./elements;
if(flag==IFFT)for(int i=0; i<elements; c[i++]*=reciprocal );
free(rotor);
}
static inline void fft_with( CZT_CONTROL *czt, double complex *c, int elements, int flag ){
if( !ft_size_check(elements) ) return; //2の冪(べき)のみ、power of two, ONLY
double complex *rotor=NULL, *table=NULL;
if( czt && czt->W[0].elements==elements && czt->W[0].flag==flag ) table=czt->W[0].table;
EI( czt && czt->W[1].elements==elements && czt->W[1].flag==flag ) table=czt->W[1].table;
if( table ) rotor = table;
EE{ double unit = (flag==IFFT) ? 2*M_PI/elements : -2*M_PI/elements;
if( !(rotor = (double complex *)calloc(elements, sizeof(double complex))) ) return;
for(int e=0; e<elements; rotor[e]=Rotor(unit*e), e++ );
}
for(int i=0, j=1; j<elements-1; j++ ) { //// scrambler, for initial placement
for(int k=elements>>1; k>(i^=k); k>>=1 );
if( j<i ){double complex temp=c[j]; c[j]=c[i], c[i]=temp; } // exchange
}
for(int m, mh=1; (m=mh<<1)<=elements; mh=m ) { //// Operations
for(int ir=0, i=0; i<elements; i+=m ) { // ir=0: Rotation angle is 0 radian.
register double complex temp, W=rotor[ir]; // W: Rotation factor
for(int k=elements>>2; k>(ir^=k); k>>=1 );
for(int j=i; j<mh+i; j++ ) { int a = j, b = j + mh;
temp = c[a] - c[b]; // For butterfly operation
c[a] = c[a] + c[b]; // a ->- (a + b) ->- a
c[b] = temp * W; // b ->- (a-b)*W ->- b
} } } double reciprocal = 1./elements;
if(flag==IFFT)for(int i=0; i<elements; c[i++]*=reciprocal );
if(!table)free(rotor);
}
////////////////////////////////////////////////////////////////////////////////
//// A function that generates a structure for Chirp Z-transform for a given
//// integer value and returns its pointer.
//// You MUST call "free()" with this pointer when it is no longer needed.
////
CZT_CONTROL *GenerateCZT_CONTROL( int n ){
int samples = n;
int samples_out = n;
int samples_fft = ft_size_check(n) ? n : get_czt_size(n);
int flag_fft = ft_size_check(n);
CZT_CONTROL *a = (CZT_CONTROL*)malloc( sizeof(CZT_CONTROL)
+ sizeof(COMPLEX)*( 2*samples_fft )
+ sizeof(COMPLEX)*( samples + (samples_fft +1)) );
// +1 for dummy
if( !a ) return NULL;
int nx = samples_fft;
*a = (CZT_CONTROL){
.samples = samples, // 標本点の数
.samples_out = samples_out, // 出力する標本点の数
.samples_fft = samples_fft, // ２の冪：２の整数乗値
.flag_fft = flag_fft, // samples が直接FFTを実行できるサイズかの判定
.W = { { nx, FFT, (double complex *)&a[1] }, // 順変換FFT用回転因子値テーブル
{ nx, IFFT, &((double complex *)&a[1])[nx] }, // 逆変換FFT用回転因子値テーブル
},
.w = &((double complex *)&a[1])[2*nx], // 重みデータ 要素数は(samples)
.v = &((double complex *)&a[1])[n+2*nx], // インパルス応答データ 要素数は(samples_fft)
};
double complex *fft_f = a->W[0].table, *fft_r = a->W[1].table;
double unit = -2*PI/nx; // for FFT
for(int e=0;e<nx;e++) fft_r[e]=conj( fft_f[e]=Rotor(unit*e) );
double complex *w=a->w, *v=a->v;
double theta = PI/n;
for(int i=n;i<samples_fft;i++) v[i]=COMPLEX(0,0);
for(int tail=samples_fft, n2=n<<1, here=0, i=0;i<n;i++ ){
double c = cos( theta*here ), s = sin( theta*here ); // 1.0*cos(), 1.0*sin() amplitude and phase
here = ( here +1 + (i<<1) ) % n2 ; // 角度情報 データー数依存なので 1例 factor.mp4
w[i] = COMPLEX( c, -s ); // complex conjugate: 0<= { w[i] <-> v[i] } < n
v[i] = COMPLEX( c, s );
v[tail-i] = COMPLEX( c, s ); // stretching ; v[tail] was keeded as dummy area
} fft_with( a, v, samples_fft, FFT ); //fft( v, samples_fft, FFT );
return a;
}
////////////////////////////////////////////////////////////////////////////////
//// Chirp Z-transform routine for fast Fourier transform of arbitrary size data
//// Argument structure is the same as the FFT routine
////
void czt( double complex *c, int n, int flag ){//For One Dimension data processing
if( ft_size_check(n) ) return fft( c, n, flag );
CZT_CONTROL *czt = NULL;
double complex *tmp = NULL; // Temporary work area used by CZT
if( !(czt = GenerateCZT_CONTROL( n ) ) ) goto GOOD_BY;
int nx = czt->samples_fft, dc_size = sizeof(double complex);
double complex *w = czt->w;
double complex *v = czt->v;
if( !(tmp = (double complex*)calloc( nx, dc_size )) ) goto GOOD_BY;
// Division is slow, so we replace it with multiplication using the reciprocal.
double m = 1.0/n; // a*m, same as a/n. if m is m = 1.0/n;
//for(int i=n;i<nx;tmp[i++]=0+I*0); // from n to nx are already zero.
if( flag == ICZT ) { // 逆変換
// ICZT part: Inverse transform processing part, inverse transform can be
// easily realized by using complex conjugate numbers.
for(int i=0;i< n;i++) tmp[i] = c[i] * conj(w[i]);
fft_with( czt, tmp, nx, FFT ); //fft( tmp, nx, FFT );
for(int i=0;i<nx;i++) tmp[i] = tmp[i] * conj(v[i]);
fft_with( czt, tmp, nx, IFFT ); //fft( tmp, nx, IFFT );
for(int i=0;i< n;i++) c[i] = tmp[i] * conj(w[i]) * m;
// Note: For the inverse transform, divide the result after the transform
// process by the number of elements.
} else {// flag==CZT 順変換
// CZT part: Forward transform processing unit
for(int i=0;i< n;i++) tmp[i] = c[i] * w[i]; // weight data multiplication
fft_with( czt, tmp, nx, FFT ); //fft( tmp, nx, FFT );
for(int i=0;i<nx;i++) tmp[i] = tmp[i] * v[i]; // convolution
fft_with( czt, tmp, nx, IFFT ); //fft( tmp, nx, IFFT );
for(int i=0;i< n;i++) c[i] = tmp[i] * w[i]; // weight data multiplication
}
GOOD_BY:;
free( tmp );
free( czt );
}