-
Notifications
You must be signed in to change notification settings - Fork 1
/
dif-org.c
115 lines (96 loc) · 2.58 KB
/
dif-org.c
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
#include "dif-org.h"
int fix_dif_fft_org(fixed fr[], fixed fi[], int m, int inverse)
{
int mr,nn,i,j,l,k,istep, n, scale, shift;
fixed qr,qi; //even input
fixed tr,ti; //odd input
fixed wr,wi; //twiddle factor
//number of input data
n = 1<<m;
if(n > N_WAVE) return -1;
mr = 0;
nn = n - 1;
scale = 0;
l = n>>1;
k = LOG2_N_WAVE-m;
while(l > 0)
{
if(inverse)
{
/* variable scaling, depending upon data */
shift = 0;
for(i=0; i<n; ++i)
{
j = fr[i];
if(j < 0) j = -j;
m = fi[i];
if(m < 0) m = -m;
if(j > 16383 || m > 16383)
{
shift = 1;
break;
}
}
if(shift) ++scale;
}
else
{
/* fixed scaling, for proper normalization -
there will be log2(n) passes, so this
results in an overall factor of 1/n,
distributed to maximize arithmetic accuracy. */
shift = 1;
}
/* it may not be obvious, but the shift will be performed
on each data point exactly once, during this pass. */
istep = l << 1; //step width of current butterfly
for(m=0; m<l; ++m)
{
j = m << k;
/* 0 <= j < N_WAVE/2 */
wr = Sinewave_org[j+N_WAVE/4];
wi = -Sinewave_org[j];
if(inverse) wi = -wi;
if(shift)
{
wr >>= 1;
wi >>= 1;
}
for(i=m; i<n; i+=istep)
{
j = i + l;
tr = fr[i] - fr[j];
ti = fi[i] - fi[j];
qr = fr[i] + fr[j];
qi = fi[i] + fi[j];
if(shift)
{
qr >>= 1;
qi >>= 1;
}
fr[i] = qr;
fi[i] = qi;
fr[j] = fix_mpy_org(wr,tr) - fix_mpy_org(wi,ti);
fi[j] = fix_mpy_org(wr,ti) + fix_mpy_org(wi,tr);
}
}
++k;
l >>= 1;
}
/* decimation in frequency - re-order data */
for(m=1; m<=nn; ++m) {
l = n;
do{
l >>= 1;
}while(mr+l > nn);
mr = (mr & (l-1)) + l;
if(mr <= m) continue;
tr = fr[m];
fr[m] = fr[mr];
fr[mr] = tr;
ti = fi[m];
fi[m] = fi[mr];
fi[mr] = ti;
}
return scale;
}