Skip to content

Commit 573841b

Browse files
agraefagraef
committed
Sanitize lyonpotpourri source file formats.
1 parent 39470d3 commit 573841b

File tree

4 files changed

+444
-4
lines changed

4 files changed

+444
-4
lines changed

‎externals/lyonpotpourri/fft.c‎

Lines changed: 160 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,160 @@
1-
#include "fftease.h"/* If forward is true, rfft replaces 2*N real data points in x with N complex values representing the positive frequency half of their Fourier spectrum, with x[1] replaced with the real part of the Nyquist frequency value. If forward is false, rfft expects x to contain a positive frequency spectrum arranged as before, and replaces it with 2*N real values. N MUST be a power of 2. */void rfft( float *x, int N, int forward ){ float c1,c2, h1r,h1i, h2r,h2i, wr,wi, wpr,wpi, temp, theta; float xr,xi; int i, i1,i2,i3,i4, N2p1; static int first = 1;/*float PI, TWOPI;*/ if ( first ) { first = 0; } theta = PI/N; wr = 1.; wi = 0.; c1 = 0.5; if ( forward ) { c2 = -0.5; cfft( x, N, forward ); xr = x[0]; xi = x[1]; } else { c2 = 0.5; theta = -theta; xr = x[1]; xi = 0.; x[1] = 0.; } wpr = -2.*pow( sin( 0.5*theta ), 2. ); wpi = sin( theta ); N2p1 = (N<<1) + 1; for ( i = 0; i <= N>>1; i++ ) { i1 = i<<1; i2 = i1 + 1; i3 = N2p1 - i2; i4 = i3 + 1; if ( i == 0 ) { h1r = c1*(x[i1] + xr ); h1i = c1*(x[i2] - xi ); h2r = -c2*(x[i2] + xi ); h2i = c2*(x[i1] - xr ); x[i1] = h1r + wr*h2r - wi*h2i; x[i2] = h1i + wr*h2i + wi*h2r; xr = h1r - wr*h2r + wi*h2i; xi = -h1i + wr*h2i + wi*h2r; } else { h1r = c1*(x[i1] + x[i3] ); h1i = c1*(x[i2] - x[i4] ); h2r = -c2*(x[i2] + x[i4] ); h2i = c2*(x[i1] - x[i3] ); x[i1] = h1r + wr*h2r - wi*h2i; x[i2] = h1i + wr*h2i + wi*h2r; x[i3] = h1r - wr*h2r + wi*h2i; x[i4] = -h1i + wr*h2i + wi*h2r; } wr = (temp = wr)*wpr - wi*wpi + wr; wi = wi*wpr + temp*wpi + wi; } if ( forward ) x[1] = xr; else cfft( x, N, forward );}/* cfft replaces float array x containing NC complex values (2*NC float values alternating real, imagininary, etc.) by its Fourier transform if forward is true, or by its inverse Fourier transform if forward is false, using a recursive Fast Fourier transform method due to Danielson and Lanczos. NC MUST be a power of 2. */void cfft( float *x, int NC, int forward ){ float wr,wi, wpr,wpi, theta, scale; int mmax, ND, m, i,j, delta; ND = NC<<1; bitreverse( x, ND ); for ( mmax = 2; mmax < ND; mmax = delta ) { delta = mmax<<1; theta = TWOPI/( forward? mmax : -mmax ); wpr = -2.*pow( sin( 0.5*theta ), 2. ); wpi = sin( theta ); wr = 1.; wi = 0.; for ( m = 0; m < mmax; m += 2 ) { register float rtemp, itemp; for ( i = m; i < ND; i += delta ) { j = i + mmax; rtemp = wr*x[j] - wi*x[j+1]; itemp = wr*x[j+1] + wi*x[j]; x[j] = x[i] - rtemp; x[j+1] = x[i+1] - itemp; x[i] += rtemp; x[i+1] += itemp; } wr = (rtemp = wr)*wpr - wi*wpi + wr; wi = wi*wpr + rtemp*wpi + wi; } }/* scale output */ scale = forward ? 1./ND : 2.; { register float *xi=x, *xe=x+ND; while ( xi < xe ) *xi++ *= scale; }}/* bitreverse places float array x containing N/2 complex values into bit-reversed order */void bitreverse( float *x, int N ){ float rtemp,itemp; int i,j, m; for ( i = j = 0; i < N; i += 2, j += m ) { if ( j > i ) { rtemp = x[j]; itemp = x[j+1]; /* complex exchange */ x[j] = x[i]; x[j+1] = x[i+1]; x[i] = rtemp; x[i+1] = itemp; } for ( m = N>>1; m >= 2 && j >= m; m >>= 1 ) j -= m; }}
1+
#include "fftease.h"
2+
3+
4+
5+
6+
/* If forward is true, rfft replaces 2*N real data points in x with
7+
N complex values representing the positive frequency half of their
8+
Fourier spectrum, with x[1] replaced with the real part of the Nyquist
9+
frequency value. If forward is false, rfft expects x to contain a
10+
positive frequency spectrum arranged as before, and replaces it with
11+
2*N real values. N MUST be a power of 2. */
12+
13+
void rfft( float *x, int N, int forward )
14+
15+
{
16+
float c1,c2,
17+
h1r,h1i,
18+
h2r,h2i,
19+
wr,wi,
20+
wpr,wpi,
21+
temp,
22+
theta;
23+
float xr,xi;
24+
int i,
25+
i1,i2,i3,i4,
26+
N2p1;
27+
static int first = 1;
28+
/*float PI, TWOPI;*/
29+
30+
if ( first ) {
31+
32+
first = 0;
33+
}
34+
theta = PI/N;
35+
wr = 1.;
36+
wi = 0.;
37+
c1 = 0.5;
38+
if ( forward ) {
39+
c2 = -0.5;
40+
cfft( x, N, forward );
41+
xr = x[0];
42+
xi = x[1];
43+
} else {
44+
c2 = 0.5;
45+
theta = -theta;
46+
xr = x[1];
47+
xi = 0.;
48+
x[1] = 0.;
49+
}
50+
wpr = -2.*pow( sin( 0.5*theta ), 2. );
51+
wpi = sin( theta );
52+
N2p1 = (N<<1) + 1;
53+
for ( i = 0; i <= N>>1; i++ ) {
54+
i1 = i<<1;
55+
i2 = i1 + 1;
56+
i3 = N2p1 - i2;
57+
i4 = i3 + 1;
58+
if ( i == 0 ) {
59+
h1r = c1*(x[i1] + xr );
60+
h1i = c1*(x[i2] - xi );
61+
h2r = -c2*(x[i2] + xi );
62+
h2i = c2*(x[i1] - xr );
63+
x[i1] = h1r + wr*h2r - wi*h2i;
64+
x[i2] = h1i + wr*h2i + wi*h2r;
65+
xr = h1r - wr*h2r + wi*h2i;
66+
xi = -h1i + wr*h2i + wi*h2r;
67+
} else {
68+
h1r = c1*(x[i1] + x[i3] );
69+
h1i = c1*(x[i2] - x[i4] );
70+
h2r = -c2*(x[i2] + x[i4] );
71+
h2i = c2*(x[i1] - x[i3] );
72+
x[i1] = h1r + wr*h2r - wi*h2i;
73+
x[i2] = h1i + wr*h2i + wi*h2r;
74+
x[i3] = h1r - wr*h2r + wi*h2i;
75+
x[i4] = -h1i + wr*h2i + wi*h2r;
76+
}
77+
wr = (temp = wr)*wpr - wi*wpi + wr;
78+
wi = wi*wpr + temp*wpi + wi;
79+
}
80+
if ( forward )
81+
x[1] = xr;
82+
else
83+
cfft( x, N, forward );
84+
}
85+
86+
/* cfft replaces float array x containing NC complex values
87+
(2*NC float values alternating real, imagininary, etc.)
88+
by its Fourier transform if forward is true, or by its
89+
inverse Fourier transform if forward is false, using a
90+
recursive Fast Fourier transform method due to Danielson
91+
and Lanczos. NC MUST be a power of 2. */
92+
93+
void cfft( float *x, int NC, int forward )
94+
95+
{
96+
float wr,wi,
97+
wpr,wpi,
98+
theta,
99+
scale;
100+
int mmax,
101+
ND,
102+
m,
103+
i,j,
104+
delta;
105+
106+
107+
ND = NC<<1;
108+
bitreverse( x, ND );
109+
for ( mmax = 2; mmax < ND; mmax = delta ) {
110+
delta = mmax<<1;
111+
theta = TWOPI/( forward? mmax : -mmax );
112+
wpr = -2.*pow( sin( 0.5*theta ), 2. );
113+
wpi = sin( theta );
114+
wr = 1.;
115+
wi = 0.;
116+
for ( m = 0; m < mmax; m += 2 ) {
117+
register float rtemp, itemp;
118+
for ( i = m; i < ND; i += delta ) {
119+
j = i + mmax;
120+
rtemp = wr*x[j] - wi*x[j+1];
121+
itemp = wr*x[j+1] + wi*x[j];
122+
x[j] = x[i] - rtemp;
123+
x[j+1] = x[i+1] - itemp;
124+
x[i] += rtemp;
125+
x[i+1] += itemp;
126+
}
127+
wr = (rtemp = wr)*wpr - wi*wpi + wr;
128+
wi = wi*wpr + rtemp*wpi + wi;
129+
}
130+
}
131+
132+
/* scale output */
133+
134+
scale = forward ? 1./ND : 2.;
135+
{ register float *xi=x, *xe=x+ND;
136+
while ( xi < xe )
137+
*xi++ *= scale;
138+
}
139+
}
140+
141+
/* bitreverse places float array x containing N/2 complex values
142+
into bit-reversed order */
143+
144+
void bitreverse( float *x, int N )
145+
146+
{
147+
float rtemp,itemp;
148+
int i,j,
149+
m;
150+
151+
for ( i = j = 0; i < N; i += 2, j += m ) {
152+
if ( j > i ) {
153+
rtemp = x[j]; itemp = x[j+1]; /* complex exchange */
154+
x[j] = x[i]; x[j+1] = x[i+1];
155+
x[i] = rtemp; x[i+1] = itemp;
156+
}
157+
for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
158+
j -= m;
159+
}
160+
}

‎externals/lyonpotpourri/leanconvert.c‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,20 @@
1-
#include "fftease.h"void leanconvert( float *S, float *C, int N2 ){ int real, imag, amp, phase; float a, b; int i; for ( i = 0; i <= N2; i++ ) { imag = phase = ( real = amp = i<<1 ) + 1; a = ( i == N2 ? S[1] : S[real] ); b = ( i == 0 || i == N2 ? 0. : S[imag] ); C[amp] = hypot( a, b ); C[phase] = -atan2( b, a ); }}
1+
#include "fftease.h"
2+
3+
void leanconvert( float *S, float *C, int N2 )
4+
5+
{
6+
7+
int real, imag,
8+
amp, phase;
9+
float a, b;
10+
int i;
11+
12+
for ( i = 0; i <= N2; i++ ) {
13+
imag = phase = ( real = amp = i<<1 ) + 1;
14+
a = ( i == N2 ? S[1] : S[real] );
15+
b = ( i == 0 || i == N2 ? 0. : S[imag] );
16+
C[amp] = hypot( a, b );
17+
C[phase] = -atan2( b, a );
18+
}
19+
}
20+

‎externals/lyonpotpourri/leanunconvert.c‎

Lines changed: 23 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,23 @@
1-
#include "fftease.h"/* unconvert essentially undoes what convert does, i.e., it turns N2+1 PAIRS of amplitude and frequency values in C into N2 PAIR of complex spectrum data (in rfft format) in output array S; sampling rate R and interpolation factor I are used to recompute phase values from frequencies */void leanunconvert( float *C, float *S, int N2 ){ int real, imag, amp, phase; register int i; for ( i = 0; i <= N2; i++ ) { imag = phase = ( real = amp = i<<1 ) + 1; S[real] = *(C+amp) * cos( *(C+phase) ); if ( i != N2 ) S[imag] = -*(C+amp) * sin( *(C+phase) ); }}
1+
#include "fftease.h"
2+
3+
/* unconvert essentially undoes what convert does, i.e., it
4+
turns N2+1 PAIRS of amplitude and frequency values in
5+
C into N2 PAIR of complex spectrum data (in rfft format)
6+
in output array S; sampling rate R and interpolation factor
7+
I are used to recompute phase values from frequencies */
8+
9+
void leanunconvert( float *C, float *S, int N2 )
10+
11+
{
12+
int real, imag,
13+
amp, phase;
14+
register int i;
15+
16+
for ( i = 0; i <= N2; i++ ) {
17+
imag = phase = ( real = amp = i<<1 ) + 1;
18+
S[real] = *(C+amp) * cos( *(C+phase) );
19+
if ( i != N2 )
20+
S[imag] = -*(C+amp) * sin( *(C+phase) );
21+
}
22+
}
23+

0 commit comments

Comments
 (0)