-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathbpf.cpp
More file actions
99 lines (85 loc) · 1.55 KB
/
bpf.cpp
File metadata and controls
99 lines (85 loc) · 1.55 KB
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
#include <math.h>
void bpf( int n, int p, float *re, float *im )
/* int n; ÞÉÓÌÏ ÔÏÞÅË - ÎÅÐÒÅÍÅÎÎÏ 2**n */
/* int p; ÔÉÐ (1 ÐÒÑÍÏÅ), (-1 ÏÂÒÁÔÎÏÅ) */
/* float re[]; ÍÁÓÓÉ× ÄÅÊÓÔ×ÉÔÅÌØÎÏÊ ÓÏÓÔÁ×ÌÑÀÝÅÊ */
/* float im[]; ÍÁÓÓÉ× ÍÎÉÍÏÊ ÓÏÓÔÁ×ÌÑÀÝÅÊ */
{
int k, l, n1, n2;
register int i, j, ii, jj;
float a, b, c, co, r, s, si, t = - 3.14159265359;
n1 = n / 2;
l = n1 / 2;
n2 = n1 + 1;
j = jj = 0;
while ( (++j) - n1 <= 0 )
{
ii = jj + n1;
r = re[j];
re[j] = re[ii];
re[ii]= r;
r = im[j];
im[j] = im[ii];
im[ii]= r;
k = l;
while ( jj - k >= 0 )
{
jj = jj - k;
k = k / 2;
}
if ( (jj += k) - (++j) <= 0 )
continue;
r = re[j];
re[j] = re[jj];
re[jj] = r;
r = im[j];
im[j] = im[jj];
im[jj] = r;
i = j + n2;
ii = jj + n2;
r = re[i];
re[i] = re[ii];
re[ii]= r;
r = im[i];
im[i] = im[ii];
im[ii]= r;
}
i = 1;
if( p < 0)
{
t = -t;
r = 1. / n;
for(j = 0; j < n; j++)
{
re[j] = re[j] * r; im[j] = im[j] * r;
}
}
do
{
si = 0.;
co = 1.;
s = sin(t);
c = cos(t);
t = 0.5 * t;
ii = i;
i = i + i;
for(k = 0; k < ii; k++)
{
for(j = k; j < n; j = j + i)
{
jj = j + ii;
a = re[jj];
b = im[jj];
r = a * co - b * si;
re[jj] = re[j] - r;
re[j] = re[j] + r;
r = b * co + a * si;
im[jj] = im[j] - r;
im[j] = im[j] + r;
}
r = c * co - s * si;
si = c * si + s * co; co = r;
}
}
while(i - n < 0);
}