File Coverage

FFT.xs
Criterion Covered Total %
statement 62 65 95.3
branch 28 32 87.5
condition n/a
subroutine n/a
pod n/a
total 90 97 92.7


line stmt bran cond sub pod time code
1             #ifdef __cplusplus
2             extern "C" {
3             #endif
4              
5             #include "EXTERN.h"
6             #include "perl.h"
7             #include "XSUB.h"
8             #include "arrays.h"
9             #include "FFT.h"
10              
11             #ifdef __cplusplus
12             }
13             #endif
14              
15             MODULE = Math::FFT PACKAGE = Math::FFT
16              
17             PROTOTYPES: DISABLE
18              
19             void
20             _cdft(n, isgn, a, ip, w)
21             int n
22             int isgn
23             double *a
24             int *ip
25             double *w
26             OUTPUT:
27             a
28              
29              
30             void
31             _rdft(n, isgn, a, ip, w)
32             int n
33             int isgn
34             double *a
35             int *ip
36             double *w
37             OUTPUT:
38             a
39              
40             void
41             _ddct(n, isgn, a, ip, w)
42             int n
43             int isgn
44             double *a
45             int *ip
46             double *w
47             OUTPUT:
48             a
49              
50             void
51             _ddst(n, isgn, a, ip, w)
52             int n
53             int isgn
54             double *a
55             int *ip
56             double *w
57             OUTPUT:
58             a
59              
60             void
61             pdfct(nt, n, a, t, ip, w)
62             int nt
63             int n
64             double *a
65             double *t
66             int *ip
67             double *w
68             CODE:
69 4           coerce1D( (SV*) ST(3), nt);
70 4           t = (double *) pack1D( (SV*) ST(3), 'd');
71 4           _dfct(n, a, t, ip, w);
72             OUTPUT:
73             a
74              
75             void
76             pdfst(nt, n, a, t, ip, w)
77             int nt
78             int n
79             double *a
80             double *t = NO_INIT
81             int *ip
82             double *w
83             CODE:
84 4           coerce1D( (SV*) ST(3), nt);
85 4           t = (double *) pack1D( (SV*) ST(3), 'd');
86 4           _dfst(n, a, t, ip, w);
87             OUTPUT:
88             a
89              
90             void
91             _correl(n, corr, d1, d2, ip, w)
92             int n
93             double *corr = NO_INIT
94             double *d1
95             double *d2
96             int *ip
97             double *w
98             PREINIT:
99             int i;
100             CODE:
101 4           coerce1D( (SV*) ST(1), n);
102 4           corr = (double *) pack1D( (SV*) ST(1), 'd');
103 4           corr[0] = d1[0]*d2[0];
104 4           corr[1] = d1[1]*d2[1];
105 32 100         for (i=2; i
106 28           corr[i] = d1[i]*d2[i]+ d1[i+1]*d2[i+1];
107 28           corr[i+1] = d1[i+1]*d2[i] - d1[i]*d2[i+1];
108             }
109 4           _rdft(n, -1, corr, ip, w);
110 68 100         for (i=0; i
111             OUTPUT:
112             corr
113              
114             void
115             _convlv(n, convlv, d1, d2, ip, w)
116             int n
117             double *convlv = NO_INIT
118             double *d1
119             double *d2
120             int *ip
121             double *w
122             PREINIT:
123             int i;
124             CODE:
125 2           coerce1D( (SV*) ST(1), n);
126 2           convlv = (double *) pack1D( (SV*) ST(1), 'd');
127 2           _rdft(n, 1, d2, ip, w);
128 2           convlv[0] = d1[0]*d2[0];
129 2           convlv[1] = d1[1]*d2[1];
130 16 100         for (i=2; i
131 14           convlv[i] = d1[i]*d2[i]- d1[i+1]*d2[i+1];
132 14           convlv[i+1] = d1[i+1]*d2[i] + d1[i]*d2[i+1];
133             }
134 2           _rdft(n, -1, convlv, ip, w);
135 34 100         for (i=0; i
136             OUTPUT:
137             convlv
138              
139              
140             int
141             _deconvlv(n, convlv, d1, d2, ip, w)
142             int n
143             double *convlv = NO_INIT
144             double *d1
145             double *d2
146             int *ip
147             double *w
148             PREINIT:
149             int i;
150             double mag;
151             CODE:
152 1           coerce1D( (SV*) ST(1), n);
153 1           convlv = (double *) pack1D( (SV*) ST(1), 'd');
154 1           _rdft(n, 1, d2, ip, w);
155 1           RETVAL=0;
156 1 50         if (fabs((float)d2[0])<1.e-10 || fabs((float)d2[1])<1.e-10) {
    50          
157 0           RETVAL=1;
158             }
159             else {
160 1           convlv[0] = d1[0]/d2[0];
161 1           convlv[1] = d1[1]/d2[1];
162 8 100         for (i=2; i
163 7           mag = d2[i]*d2[i] + d2[i+1]*d2[i+1];
164 7 50         if (mag < 1.0e-10) {
165 0           RETVAL =1;
166 0           break;
167             }
168 7           convlv[i] = (d1[i]*d2[i]+ d1[i+1]*d2[i+1])/mag;
169 7           convlv[i+1] = (d1[i+1]*d2[i] - d1[i]*d2[i+1])/mag;
170             }
171 1 50         if (RETVAL == 0) {
172 1           _rdft(n, -1, convlv, ip, w);
173 17 100         for (i=0; i
174             }
175             }
176             OUTPUT:
177             convlv
178             RETVAL
179              
180             void
181             _spctrm(n, spctrm, data, ip, w, n2, flag)
182             int n
183             double *spctrm = NO_INIT
184             double *data
185             int *ip
186             double *w
187             double n2
188             int flag
189             PREINIT:
190             int i;
191             CODE:
192 2           coerce1D( (SV*) ST(1), n/2+1);
193 2           spctrm = (double *) pack1D( (SV*) ST(1), 'd');
194 2 100         if (flag == 1) _rdft(n, 1, data, ip, w);
195 2           spctrm[0] = data[0]*data[0] / n2;
196 2           spctrm[n/2] = data[1]*data[1] / n2;
197 16 100         for (i=1; i
198 14           spctrm[i] = 2*(data[2*i]*data[2*i]+ data[2*i+1]*data[2*i+1])/n2;
199             }
200             OUTPUT:
201             spctrm
202              
203             void
204             _spctrm_bin(k, m, spctrm, data, ip, w, n2, tmp)
205             int k
206             int m
207             double *spctrm = NO_INIT
208             double2D *data
209             int *ip
210             double *w
211             double n2
212             double *tmp = NO_INIT
213             PREINIT:
214             int i, j, n;
215 10           double den = 0.0;
216             CODE:
217 10           coerce1D( (SV*) ST(2), m/2+1);
218 10           spctrm = (double *) pack1D( (SV*) ST(2), 'd');
219 10           coerce1D( (SV*) ST(7), m);
220 10           tmp = (double *) pack1D( (SV*) ST(7), 'd');
221 250 100         for(i=0; i
222 7920 100         for (j=0; j
223 240           _rdft(m, 1, tmp, ip, w);
224 240           spctrm[0] += tmp[0]*tmp[0];
225 240           spctrm[m/2] += tmp[1]*tmp[1];
226 240           den += n2;
227 3840 100         for (n=1; n
228 3600           spctrm[n] += 2*(tmp[2*n]*tmp[2*n]+ tmp[2*n+1]*tmp[2*n+1]);
229             }
230 10           den *= m;
231 180 100         for (j=0; j<=m/2; j++) spctrm[j] /= den;
232             OUTPUT:
233             spctrm
234