#ifdef __cplusplus
extern "C" {
#endif
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "arrays.h"
#include "FFT.h"
#ifdef __cplusplus
}
#endif
MODULE = Math::FFT PACKAGE = Math::FFT
PROTOTYPES: DISABLE
void
_cdft(n, isgn, a, ip, w)
int n
int isgn
double *a
int *ip
double *w
OUTPUT:
a
void
_rdft(n, isgn, a, ip, w)
int n
int isgn
double *a
int *ip
double *w
OUTPUT:
a
void
_ddct(n, isgn, a, ip, w)
int n
int isgn
double *a
int *ip
double *w
OUTPUT:
a
void
_ddst(n, isgn, a, ip, w)
int n
int isgn
double *a
int *ip
double *w
OUTPUT:
a
void
pdfct(nt, n, a, t, ip, w)
int nt
int n
double *a
double *t
int *ip
double *w
CODE:
coerce1D( (SV*) ST(3), nt);
t = (double *) pack1D( (SV*) ST(3), 'd');
_dfct(n, a, t, ip, w);
OUTPUT:
a
void
pdfst(nt, n, a, t, ip, w)
int nt
int n
double *a
double *t = NO_INIT
int *ip
double *w
CODE:
coerce1D( (SV*) ST(3), nt);
t = (double *) pack1D( (SV*) ST(3), 'd');
_dfst(n, a, t, ip, w);
OUTPUT:
a
void
_correl(n, corr, d1, d2, ip, w)
int n
double *corr = NO_INIT
double *d1
double *d2
int *ip
double *w
PREINIT:
int i;
CODE:
coerce1D( (SV*) ST(1), n);
corr = (double *) pack1D( (SV*) ST(1), 'd');
corr[0] = d1[0]*d2[0];
corr[1] = d1[1]*d2[1];
for (i=2; i<n; i+=2) {
corr[i] = d1[i]*d2[i]+ d1[i+1]*d2[i+1];
corr[i+1] = d1[i+1]*d2[i] - d1[i]*d2[i+1];
}
_rdft(n, -1, corr, ip, w);
for (i=0; i<n; i++) corr[i] *= 2.0/n;
OUTPUT:
corr
void
_convlv(n, convlv, d1, d2, ip, w)
int n
double *convlv = NO_INIT
double *d1
double *d2
int *ip
double *w
PREINIT:
int i;
CODE:
coerce1D( (SV*) ST(1), n);
convlv = (double *) pack1D( (SV*) ST(1), 'd');
_rdft(n, 1, d2, ip, w);
convlv[0] = d1[0]*d2[0];
convlv[1] = d1[1]*d2[1];
for (i=2; i<n; i+=2) {
convlv[i] = d1[i]*d2[i]- d1[i+1]*d2[i+1];
convlv[i+1] = d1[i+1]*d2[i] + d1[i]*d2[i+1];
}
_rdft(n, -1, convlv, ip, w);
for (i=0; i<n; i++) convlv[i] *= 2.0/n;
OUTPUT:
convlv
int
_deconvlv(n, convlv, d1, d2, ip, w)
int n
double *convlv = NO_INIT
double *d1
double *d2
int *ip
double *w
PREINIT:
int i;
double mag;
CODE:
coerce1D( (SV*) ST(1), n);
convlv = (double *) pack1D( (SV*) ST(1), 'd');
_rdft(n, 1, d2, ip, w);
RETVAL=0;
if (fabs((float)d2[0])<1.e-10 || fabs((float)d2[1])<1.e-10) {
RETVAL=1;
}
else {
convlv[0] = d1[0]/d2[0];
convlv[1] = d1[1]/d2[1];
for (i=2; i<n; i+=2) {
mag = d2[i]*d2[i] + d2[i+1]*d2[i+1];
if (mag < 1.0e-10) {
RETVAL =1;
break;
}
convlv[i] = (d1[i]*d2[i]+ d1[i+1]*d2[i+1])/mag;
convlv[i+1] = (d1[i+1]*d2[i] - d1[i]*d2[i+1])/mag;
}
if (RETVAL == 0) {
_rdft(n, -1, convlv, ip, w);
for (i=0; i<n; i++) convlv[i] *= 2.0/n;
}
}
OUTPUT:
convlv
RETVAL
void
_spctrm(n, spctrm, data, ip, w, n2, flag)
int n
double *spctrm = NO_INIT
double *data
int *ip
double *w
int n2
int flag
PREINIT:
int i;
CODE:
coerce1D( (SV*) ST(1), n/2+1);
spctrm = (double *) pack1D( (SV*) ST(1), 'd');
if (flag == 1) _rdft(n, 1, data, ip, w);
spctrm[0] = data[0]*data[0] / n2;
spctrm[n/2] = data[1]*data[1] / n2;
for (i=1; i<n/2; i++) {
spctrm[i] = 2*(data[2*i]*data[2*i]+ data[2*i+1]*data[2*i+1])/n2;
}
OUTPUT:
spctrm
void
_spctrm_bin(k, m, spctrm, data, ip, w, n2, tmp)
int k
int m
double *spctrm = NO_INIT
double2D *data
int *ip
double *w
double n2
double *tmp = NO_INIT
PREINIT:
int i, j, n;
double den = 0.0;
CODE:
coerce1D( (SV*) ST(2), m/2+1);
spctrm = (double *) pack1D( (SV*) ST(2), 'd');
coerce1D( (SV*) ST(7), m);
tmp = (double *) pack1D( (SV*) ST(7), 'd');
for(i=0; i<k*m; i+=m) {
for (j=0; j<m; j++) tmp[j] = data[i+j];
_rdft(m, 1, tmp, ip, w);
spctrm[0] += tmp[0]*tmp[0];
spctrm[m/2] += tmp[1]*tmp[1];
den += n2;
for (n=1; n<m/2; n++)
spctrm[n] += 2*(tmp[2*n]*tmp[2*n]+ tmp[2*n+1]*tmp[2*n+1]);
}
den *= m;
for (j=0; j<=m/2; j++) spctrm[j] /= den;
OUTPUT:
spctrm