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
|
/* Date: Fri, 26 Sep 1997 23:26:27 +0200
From: Julia Lawall <Julia.Lawall@irisa.fr>
Subject: FFT program
The following is a very naive implementation that most of my paper
centers around. I got it from the following:
@Misc{SimplFFT,
author = {Arndt, J.},
howpublished = {{\tt
http://www.spektracom.de/\~{}arndt/fxt/fxt970326.tgz} in
the file {\tt hfloat/src/fxt/simplfft/fft.c}}
}
*/
/* from the ~/spec/fxt/hfloat/hfloat/src/fxt/simplfft/fft.c
with some simplifications to get rid of the realfft implementation */
/* Has been adjusted to C-Mix/II by Arne John Glenstrup <panic@diku.dk> */
#define union struct
#include <math.h>
#include "fft.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define SWAP(x,y) {REAL tmp=x; x=y; y=tmp;}
extern REAL orig_real[N], orig_imag[N];
REAL real[N], imag[N];
void scramble(int n2, REAL *real, REAL *imag)
{
long m,j;
for (m=1,j=0; m<n2-1; m++)
{
long k;
{ for(k=n2>>1; (!((j^=k)&k));) k>>=1; }
if (j>m)
{
SWAP(real[m],real[j]);
SWAP(imag[m],imag[j]);
}
}
}
void fft(int is, REAL *fr, REAL *fi)
{
int n2;
int ldm,m,mh,j;
int r;
int t1, t2;
REAL pi,phi,c,s,w;
REAL ur,vr, ui,vi;
n2 = 1<<LOGN;
pi = is * M_PI;
scramble(n2, fr, fi);
for (ldm=1; ldm<=LOGN; ldm++)
{
m = (1<<ldm);
mh = (m>>1);
phi = pi/(REAL)(mh);
for (j=0; j<mh; j++)
{
w = phi * (REAL)j;
c = cos(w);
s = sin(w);
for (r=0; r<n2; )
{
t1 = r + j;
t2 = t1 + mh;
vr = fr[t2] * c - fi[t2] * s;
vi = fr[t2] * s + fi[t2] * c;
ur = fr[t1];
fr[t1] = fr[t1] + vr;
fr[t2] = ur - vr;
ui = fi[t1];
fi[t1] = fi[t1] + vi;
fi[t2] = ui - vi;
r = r + m;
}
}
}
}
/* ---------------------------------------------------------------------- */
/* Apply the fft to the generated test data, based on q, iter times */
REAL applyfft() {
int i;
int x;
REAL total;
total = 0.0;
/* reset the real and imag arrays to the generated test data */
for (x = 0; x < N; x++) {
real[x] = orig_real[x];
imag[x] = orig_imag[x];
}
fft(1, real, imag);
for (x = 0; x < N; x++) {
total = total + real[x];
total = total + imag[x];
}
return total;
}
|
