-
Notifications
You must be signed in to change notification settings - Fork 0
/
FFT.xs
234 lines (216 loc) · 4.2 KB
/
FFT.xs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
#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