Two dimentional DFT

twoDimDFTTest.cpp

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
using namespace std;

typedef struct _complex {
  double re;
  double im;
} complex;


void **calloc2(size_t size, int x, int y);
void free2(int size, void **tgt);
void initializeTrigonometricTable(int size, double **tcos, double **tsin);
void oneDimDFT(int n,
	       complex *base, complex *dft,
	       double **tcos, double **tsin
	       );
void oneDimIDFT(int n,
		complex *dft, complex *idft,
		double **tcos, double **tsin
		);
void twoDimDFT(int x, int y,
	       complex **base, complex **dft,
	       double **tcosX, double **tsinX,
	       double **tcosY, double **tsinY
	       );
void twoDimIDFT(int x, int y,
	       complex **dft, complex **idft,
	       double **tcosX, double **tsinX,
	       double **tcosY, double **tsinY
	       );
void printComplex1(int n, complex *tgt);
void printComplex2(int x, int y, complex **table);

int main(void){
  int x = 2;
  int y = 3;
  double
    **tcosX = NULL, **tsinX = NULL,
    **tcosY = NULL, **tsinY = NULL;
  tcosX = (double **)calloc2(sizeof(double), x, x);
  tsinX = (double **)calloc2(sizeof(double), x, x);
  assert(tcosX != NULL && tsinX != NULL);
  initializeTrigonometricTable(x, tcosX, tsinX);
  tcosY = (double **)calloc2(sizeof(double), y, y);
  tsinY = (double **)calloc2(sizeof(double), y, y);
  assert(tcosY != NULL && tsinY != NULL);
  initializeTrigonometricTable(y, tcosY, tsinY);


  complex *oneDim, *odft, *oidft;
  cout << "==One dimensional dft & idft test==" << endl;
  oneDim = (complex*)malloc(x*sizeof(complex));
  odft = (complex*)malloc(x*sizeof(complex));
  oidft = (complex*)malloc(x*sizeof(complex));
  assert(oneDim != NULL && odft != NULL && oidft != NULL);
  for(int i = 0; i < x; i++){
    oneDim[i].re = i+1;
    oneDim[i].im = (x-i)+((double)(y-i)/10.0);
  }
  printComplex1(x, oneDim);
  oneDimDFT(x, oneDim, odft, tcosX, tsinX);
  oneDimIDFT(x, odft, oidft, tcosX, tsinX);
  printComplex1(x, oidft);
  free(oneDim);free(odft);free(oidft);
  cout << "==end==" << endl << endl;
  
  complex **base, **dft , **idft;
  base = (complex **)calloc2(sizeof(complex), y, x);
  dft = (complex **)calloc2(sizeof(complex), y, x);
  idft = (complex **)calloc2(sizeof(complex), y, x);
  
  int cnt = 1;
  for (int j = 0; j < y; j++) {
    for (int i = 0; i < x; i++) {
      base[j][i].re = (double)cnt++;
      base[j][i].im = (double)(x*y)-cnt;
    }
  }

  cout << "==Two dimensional dft & idft test==" << endl;
  cout << "target" << endl;
  printComplex2(x, y, base);
  twoDimDFT(x, y, base, dft, tcosX, tsinX, tcosY, tsinY);
  cout << endl << "dft" << endl;
  printComplex2(x, y, dft);
  twoDimIDFT(x, y, dft, idft, tcosX, tsinX, tcosY, tsinY);
  cout << endl << "idft" << endl;
  printComplex2(x, y, idft);
  cout << "==end==" << endl;

  free2(x, (void **)tcosX);free2(x, (void **)tsinX);
  free2(y, (void **)tcosY);free2(y, (void **)tsinY);
  free2(y, (void **)base);
  free2(y, (void **)dft);
  free2(y, (void **)idft);
  return 0;
}

void **calloc2(size_t size, int x, int y)
{
  void **arr = (void **)malloc(x * sizeof(void *));
  assert(arr != NULL);
  for (int i = 0; i < x; i++) {
    arr[i] = calloc(y, size);
  }
  return arr;
}

void free2(int size, void **tgt)
{
  for (int i = 0; i < size; i++) {
    free(tgt[i]);
  }
  free(tgt);
}

void initializeTrigonometricTable(int size, double **tcos, double **tsin)
{
  assert(size > 1 && tcos != NULL && *tsin != NULL);
  for (int k = 0; k < size; k++ ) {
    for (int i = 0; i < size; i++) {
      double radian = 2.0 * M_PI / (double)size * (double)k * (double)i;
      tcos[k][i] = cos(radian);
      tsin[k][i] = sin(radian);
    }
  }
}

void printComplex1(int n, complex *tgt)
{
  printf("%3s %10s\n", "num", "value");
  for (int i = 0; i < n; i++) {
    printf("%3d  %2.3lf%+2.3lfi\n", i, tgt[i].re, tgt[i].im);
  }
}

void printComplex2(int x, int y, complex **table)
{
  for (int j = 0; j<y; j++ ) {
    for (int i = 0; i < x; i++) {
      printf("%2.3lf%+2.3lfi ", table[j][i].re, table[j][i].im);
    }
    cout << endl;
  }
}

void oneDimDFT(int n,
	       complex *base, complex *dft,
	       double **tcos, double **tsin
	       )
{
  double sum_of_re;
  double sum_of_im;
  
  for (int k = 0; k < n; k++) {
    sum_of_re = 0;
    sum_of_im = 0;    
    for (int i = 0; i < n; i++) {
      sum_of_re += base[i].re * tcos[k][i] + base[i].im * tsin[k][i];
      sum_of_im += base[i].im * tcos[k][i] - base[i].re * tsin[k][i];
    }
    dft[k].re = sum_of_re / n;
    dft[k].im = sum_of_im / n;
  }
}

void oneDimIDFT(int n,
		complex *dft, complex *idft,
		double **tcos, double **tsin
		)
{
  double sum_of_re;
  double sum_of_im;

  for (int i = 0; i < n; i++) {
    sum_of_re = 0;
    sum_of_im = 0;
    for (int k = 0; k < n; k++) {
      sum_of_re += dft[k].re * tcos[k][i] - dft[k].im * tsin[k][i];
      sum_of_im += dft[k].im * tcos[k][i] + dft[k].re * tsin[k][i];
    }
    idft[i].re = sum_of_re;
    idft[i].im = sum_of_im;
  }
}

void twoDimDFT(int x, int y,
	       complex **base, complex **dft,
	       double **tcosX, double **tsinX,
	       double **tcosY, double **tsinY
	       )
{
  // base, dft must $(name)[y][x]
  complex *tmp1, *tmp2;

  // X
  tmp1 = (complex*)calloc(x, sizeof(complex));
  tmp2 = (complex*)calloc(x, sizeof(complex));
  for (int j = 0; j < y; j++) {
    for (int i = 0; i < x; i++) {
      tmp1[i] = base[j][i];
    }
    oneDimDFT(x, tmp1, tmp2, tcosX, tsinX);
    for (int i = 0; i < x; i++) {
      dft[j][i] = tmp2[i];
    }
  }
  free(tmp1); free(tmp2);


  // Y
  tmp1 = (complex*)calloc(y, sizeof(complex));
  tmp2 = (complex*)calloc(y, sizeof(complex));
  for (int i = 0; i < x; i++) {
    for (int j = 0; j < y; j++) {
      tmp1[j] = dft[j][i];
    }
    oneDimDFT(y, tmp1, tmp2, tcosY, tsinY);
    for (int j = 0; j < y; j++) {
      dft[j][i] = tmp2[j];
    }
  }
  free(tmp1); free(tmp2);
}

void twoDimIDFT(int x, int y,
	       complex **dft, complex **idft,
	       double **tcosX, double **tsinX,
	       double **tcosY, double **tsinY
	       )
{
  // dft, idft must $(name)[y][x]
  complex *tmp1, *tmp2;
  
  // Y
  tmp1 = (complex*)calloc(y, sizeof(complex));
  tmp2 = (complex*)calloc(y, sizeof(complex));
  for (int i = 0; i < x; i++) {
    for (int j = 0; j < y; j++) {
      tmp1[j] = dft[j][i];
    }
    oneDimIDFT(y, tmp1, tmp2, tcosY, tsinY);
    for (int j = 0; j < y; j++) {
      idft[j][i] = tmp2[j];
    }
  }
  free(tmp1); free(tmp2);

  // X
  tmp1 = (complex*)calloc(x, sizeof(complex));
  tmp2 = (complex*)calloc(x, sizeof(complex));
  for (int j = 0; j < y; j++) {
    for (int i = 0; i < x; i++) {
      tmp1[i] = idft[j][i];
    }
    oneDimIDFT(x, tmp1, tmp2, tcosX, tsinX);
    for (int i = 0; i < x; i++) {
      idft[j][i] = tmp2[i];
    }
  }
  free(tmp1); free(tmp2);
}