From 5b7884bf5ac56324ff52a3e427ef68c350193b09 Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Wed, 10 Jun 2020 19:49:52 +0200
Subject: [PATCH] fftlib
---
FFTlib/wkykx2yxt.c | 113 +++++++++++++++++++++++++++++++++++++++++++++
FFTlib/yxt2wkykx.c | 108 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 221 insertions(+)
create mode 100644 FFTlib/wkykx2yxt.c
create mode 100644 FFTlib/yxt2wkykx.c
diff --git a/FFTlib/wkykx2yxt.c b/FFTlib/wkykx2yxt.c
new file mode 100644
index 0000000..109a4aa
--- /dev/null
+++ b/FFTlib/wkykx2yxt.c
@@ -0,0 +1,113 @@
+#include "genfft.h"
+
+/**
+NAME: wkx2xt
+
+DESCRIPTION:
+ 2 Dimensional complex to real FFT (omega,kx -> x,t) with scaling
+
+USAGE:
+ void wkx2xt(complex *cdata, REAL *rdata, int nt, int nx,
+ int ldc, int ldr, int xorig)
+
+INPUT:
+ - *cdata: complex 2D input array [nf][nkx]
+ - nt: number of real (fast) output samples
+ - nx: number of complex (slow) samples to be transformed
+ - ldc: leading dimension (number of complex samples)
+ - ldr: leading dimension (number of real samples)
+ - xorig: trace number of origin of x-axis first trace # 0
+
+OUTPUT: - *rdata: real 2D output array scaled [nx][nt]
+
+AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+----------------------------------------------------------------------*/
+
+void wkykx2yxt(complex *cdata, REAL *rdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig)
+{
+ int i, j, k, ld1, sign, nxorig, nyorig;
+ REAL scl;
+ complex *cdum;
+
+
+ assert ( (xorig >= 0) && (xorig < nx) );
+ assert ( (yorig >= 0) && (yorig < ny) );
+ scl = 1.0/(nt*nx*ny);
+ ld1 = (nt+2)/2;
+ nxorig = nx-xorig;
+ nyorig = ny-yorig;
+
+ cdum = (complex *)malloc(ld1*ldx*ldy*sizeof(complex));
+
+ sign = -1;
+ // ctrans [nkx][nf][nky]
+ for (i = 0; i < ldx; i++) {
+ for (j = 0; j < ld1; j++) {
+ for (k = 0; k < ldy; k++) {
+ cdum[i*ld1*ldy+j*ldy+k] = cdata[j*ldy*ldx+k*ldx+i];
+ }
+ }
+ ccmfft(&cdum[i*ld1*ldy], ny, ld1, ldy, sign);
+ }
+
+ // cdata [nf][nky][nkx]
+ for (k = 0; k < ldy; k++) {
+ for (j = 0; j < ld1; j++) {
+ for (i = 0; i < ldx; i++) {
+ cdata[k*ld1*ldx+j*ldx+i] = cdum[i*ld1*ldy+j*ldy+k];
+ }
+ }
+ ccmfft(&cdata[k*ld1*ldx], nx, ld1, ldx, sign);
+ }
+ free(cdum);
+
+ cdum = (complex *)malloc(ld1*nx*ny*sizeof(complex));
+ if (cdum == NULL) fprintf(stderr,"wkx2xt: memory allocation error\n");
+
+ //cdum [ny][nx][nf]
+ for (j = 0; j < ld1; j++) {
+ for (k = 0; k < nyorig; k++) {
+ for (i = 0; i < nxorig; i++) {
+ cdum[(k+yorig)*nx*ld1+(i+xorig)*ld1+j] = cdata[k*ld1*ldx+j*ldx+i];
+ }
+ for(i = 0; i < xorig; i++) {
+ cdum[(k+yorig)*nx*ld1+i*ld1+j] = cdata[k*ld1*ldx+j*ldx+i+nxorig];
+ }
+ }
+ for (k = 0; k < yorig; k++) {
+ for (i = 0; i < nxorig; i++) {
+ cdum[k*nx*ld1+(i+xorig)*ld1+j] = cdata[k*ld1*ldx+j*ldx+i];
+ }
+ for(i = 0; i < xorig; i++) {
+ cdum[k*nx*ld1+i*ld1+j] = cdata[k*ld1*ldx+j*ldx+i+nxorig];
+ }
+ }
+ }
+
+ sign = 1;
+ crmfft(cdum, rdata, nt, nx*ny, ld1, ldt, sign);
+
+ for(i = 0; i < nx*ny*ldt; i++) rdata[i] *= scl;
+
+ free(cdum);
+
+ return;
+}
+
+/****************** FORTRAN SHELL *****************/
+
+#ifdef DF_CAPFNAMES
+#define nwkykx2yxt FNAME(WKYKX2YXTF)
+#else
+#define nwkykx2yxt FNAME(wkykx2yxtf)
+#endif
+
+void nwkykx2yxt(complex *cdata, REAL *rdata, int *nt, int *nx, int *ny, int *ldt, int *ldx, int *ldy, int *xorig, int *yorig)
+{
+
+ wkykx2yxt(cdata, rdata, *nt, *nx, *ny, *ldt, *ldx, *ldy, *xorig, *yorig);
+
+ return;
+}
diff --git a/FFTlib/yxt2wkykx.c b/FFTlib/yxt2wkykx.c
new file mode 100644
index 0000000..53a236b
--- /dev/null
+++ b/FFTlib/yxt2wkykx.c
@@ -0,0 +1,108 @@
+#include "genfft.h"
+
+/**
+NAME: xt2wkx
+
+DESCRIPTION:
+ 2 Dimensional real to complex FFT (x,t -> omega,kx) with scaling
+
+USAGE:
+ void xt2wkx(REAL *rdata, complex *cdata, int nt, int nx,
+ int ldr, int ldc, int xorig)
+
+INPUT:
+ - *rdata: real 2D input array [ny][nx][nt]
+ - nt: number of real (fast) output samples
+ - nx: number of complex (slow) samples to be transformed
+ - ny: number of complex (slow) samples to be transformed
+ - ldt: leading dimension (number of real samples)
+ - ldx: leading dimension (number of complex samples)
+ - ldy: leading dimension (number of complex samples)
+ - xorig: trace number of origin of x-axis first trace # 0
+ - yorig: trace number of origin of y-axis first trace # 0
+
+OUTPUT: - *cdata: complex 2D output array scaled [nf][nky][nkx]
+
+AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+----------------------------------------------------------------------*/
+
+void yxt2wkykx(REAL *rdata, complex *cdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig)
+{
+ long i, j, k, ld1, sign;
+ complex *cdum;
+
+ assert ( (xorig >= 0) && (xorig < nx) );
+ assert ( (yorig >= 0) && (yorig < ny) );
+ ld1 = (nt+2)/2;
+ cdum = (complex *)malloc(ld1*nx*ny*sizeof(complex));
+ if (cdum == NULL) fprintf(stderr,"yxt2wkykx: memory allocation error\n");
+
+ sign = -1;
+ rcmfft(rdata, cdum, nt, nx*ny, ldt, ld1, sign);
+
+ // cdata[nky][nf][nkx] = cdum[ny][nx][nf]
+ for(j = 0; j < ld1; j++) {
+ for(k = yorig; k < ny; k++) {
+ for(i = xorig; i < nx; i++) {
+ cdata[(k-yorig)*ld1*ldx+j*ldx+(i-xorig)] = cdum[k*nx*ld1+i*ld1+j];
+ }
+ for(i = 0; i < xorig; i++) {
+ cdata[(k-yorig)*ld1*ldx+j*ldx+(nx-xorig+i)] = cdum[k*nx*ld1+i*ld1+j];
+ }
+ }
+ for(k = 0; k < yorig; k++) {
+ for(i = xorig; i < nx; i++) {
+ cdata[(ny-yorig+k)*ld1*ldx+j*ldx+(i-xorig)] = cdum[k*nx*ld1+i*ld1+j];
+ }
+ for(i = 0; i < xorig; i++) {
+ cdata[(ny-yorig+k)*ld1*ldx+j*ldx+(nx-xorig+i)] = cdum[k*nx*ld1+i*ld1+j];
+ }
+ }
+ }
+
+ free(cdum);
+
+ cdum = (complex *)malloc(ld1*ldx*ldy*sizeof(complex));
+
+ sign = 1;
+ // cdum[nkx][nf][nky] = cdata[nky][nf][nkx]
+ for (k = 0; k < ldy; k++) {
+ ccmfft(&cdata[k*ld1*ldx], nx, ld1, ldx, sign);
+ for (j = 0; j < ld1; j++) {
+ for (i = 0; i < ldx; i++) {
+ cdum[i*ld1*ldy+j*ldy+k] = cdata[k*ld1*ldx+j*ldx+i];
+ }
+ }
+ }
+
+ // cdata [nf][nky][nkx] = cdum[nkx][nf][nky]
+ for (i = 0; i < ldx; i++) {
+ ccmfft(&cdum[i*ld1*ldy], ny, ld1, ldy, sign);
+ for (j = 0; j < ld1; j++) {
+ for (k = 0; k < ldy; k++) {
+ cdata[j*ldy*ldx+k*ldx+i] = cdum[i*ld1*ldy+j*ldy+k];
+ }
+ }
+ }
+ free(cdum);
+
+ return;
+}
+
+/****************** FORTRAN SHELL *****************/
+
+#ifdef DF_CAPFNAMES
+#define nynxt2wkykx FNAME(YXT2WKYKXF)
+#else
+#define nynxt2wkykx FNAME(yxt2wkykxf)
+#endif
+
+void nynxt2wkykx(REAL *rdata, complex *cdata, int *nt, int *nx, int *ny, int *ldt, int *ldx, int *ldy, int *xorig, int *yorig)
+{
+
+ yxt2wkykx(rdata, cdata, *nt, *nx, *ny, *ldt, *ldx, *ldy, *xorig, *yorig);
+
+ return;
+}
--
GitLab