1 // This file is part of the FXT library.
2 // Copyright (C) 2010, 2011, 2012 Joerg Arndt
3 // License: GNU General Public License version 3 or later,
4 // see the file COPYING.txt in the main directory.
5
6 #include "chirpzt/chirpzt.h"
7 #include "convolution/complexconvolution.h"
8
9 #include "bits/bit2pow.h" // ld()
10 #include "aux1/arith1.h" // negate()
11
12 //#include "aux1/aux1double.h" // ri_multiply()
13 #include "aux0/cmult.h" // cmult(), used in ri_multiply()
14
15 #include "aux0/sincos.h"
16 #include "aux0/csincos.h"
17
18 #include "aux1/copy.h"
19
20 #include "fxttypes.h" // ulong
21 #include "complextype.h" // Complex
22
23 #include <cmath> // M_PI
24
25
26
27 static inline void
make_fft_fract_chirp(Complex * w,double v,ulong n,ulong nn)28 make_fft_fract_chirp(Complex *w, double v, ulong n, ulong nn)
29 // For k=0..nn-1: w[k] == exp(v*sqrt(-1)*k*k*2*pi*/n/2)
30 {
31 // const double phi = v*2.0*M_PI/(double)n/2.0;
32 const double phi = v*M_PI/(double)n;
33 ulong n2 = 2*n;
34 ulong np = 0;
35 for (ulong k=0; k<nn; ++k)
36 {
37 w[k] = SinCos(phi*(double)np);
38
39 np += ((k<<1)+1); // np == (k*k)%n2
40 if ( np>=n2 ) np -= n2;
41 }
42 }
43 // -------------------------
44
45
46 void
fft_fract(Complex * x,ulong n,double v)47 fft_fract(Complex *x, ulong n, double v)
48 // Fractional (fast) Fourier transform.
49 //
50 // For complex array c[0...n]
51 // compute \sum_{x=0}^{n}{c_x*exp(is*v*2*i*\pi*x*k/n)}
52 // (for v==1.0 this is just the usual FFT)
53 //
54 // Use n*k == n^2/2 + k^2/2 - (k-n)^2/2
55 // See: Nussbaumer: "Fast Fourier Transform and Convolution Algorithms", chap.5, p.113ff
56 // Could use n*k == - n^2/2 - k^2/2 + (k+n)^2/2 instead.
57 //
58 // Worst case if n=2^x+1:
59 // then nn=2*2^x
60 // work is about 6 times of an FFT of length 2^x
61 // and allocated workspace =2*nn
62 {
63 const ulong ldnn = 1 + ld((n<<1)-1);
64 const ulong nn = (1UL<<ldnn); // smallest power of 2 >= 2*n
65
66 Complex *f = new Complex[nn];
67 acopy(x, f, n);
68 null(f+n, nn-n);
69
70 Complex *w = new Complex[nn];
71 make_fft_fract_chirp(w, v, n, nn);
72
73 // multiply(w, f, n);
74 for (ulong j=0; j<n; ++j) f[j] *= w[j];
75
76 // negate(wi, nn);
77 for (ulong j=0; j<nn; ++j) w[j] = conj(w[j]);
78
79 fft_complex_convolution(w, f, ldnn);
80
81 make_fft_fract_chirp(w, v, n, nn);
82
83 // multiply(f, w, nn);
84 for (ulong j=0; j<n; ++j) w[j] *= f[j];
85
86 acopy(w+n, x, n);
87 delete [] w;
88 delete [] f;
89 }
90 // -------------------------
91
92
93
94 static inline void
make_fft_fract_chirp(double * wr,double * wi,double v,ulong n,ulong nn)95 make_fft_fract_chirp(double *wr, double *wi, double v, ulong n, ulong nn)
96 // For k=0..nn-1: Complex(wr[k], wi[k]) == exp(v*sqrt(-1)*k*k*2*pi*/n/2)
97 {
98 // const double phi = v*2.0*M_PI/n/2;
99 const double phi = v*M_PI/(double)n;
100 ulong n2 = 2*n;
101 ulong np=0;
102 for (ulong k=0; k<nn; ++k)
103 {
104 double c, s;
105 SinCos(phi*(double)np, &s, &c);
106 wr[k] = c;
107 wi[k] = s;
108
109 np += ((k<<1)+1); // np == (k*k)%n2
110 if ( np>=n2 ) np -= n2;
111 }
112 }
113 // -------------------------
114
115
116
ri_multiply(const double * fr,const double * fi,double * gr,double * gi,ulong n)117 static void ri_multiply(const double *fr, const double *fi,
118 double *gr, double *gi, ulong n)
119 // complex multiply array g[] by f[]:
120 // Complex(gr[], gi[]) *= Complex(fr[], fi[])
121 {
122 while ( n-- ) cmult(fr[n], fi[n], gr[n], gi[n]);
123 }
124 // -------------------------
125
126 void
fft_fract(double * x,double * y,ulong n,double v)127 fft_fract(double *x, double *y, ulong n, double v)
128 // Fractional (fast) fourier transform.
129 //
130 // For complex array c[0...n]
131 // compute \sum_{x=0}^{n}{c_x*exp(is*v*2*i*\pi*x*k/n)}
132 // (for v==1.0 this is just the usual fft)
133 //
134 // Use n*k == n^2/2 + k^2/2 - (k-n)^2/2
135 // See: Nussbaumer: "Fast Fourier Transform and Convolution Algorithms", chap.5, p.113ff
136 // Could use n*k == - n^2/2 - k^2/2 + (k+n)^2/2 instead.
137 //
138 // Worst case if n=2^x+1:
139 // then nn=4*2^x
140 // work is about 12 times of an FFT of length 2^x
141 // and allocated workspace =4*nn
142 {
143 const ulong ldnn = 1 + ld((n<<1)-1);
144 const ulong nn = (1UL<<ldnn); // smallest power of 2 >= 2*n
145
146 double *fr = new double[2*nn];
147 double *fi = fr + nn;
148
149 acopy(x, fr, n);
150 acopy(y, fi, n);
151
152 null(fr+n, nn-n);
153 null(fi+n, nn-n);
154
155 double *wr = new double[2*nn];
156 double *wi = wr + nn;
157
158 make_fft_fract_chirp(wr, wi, v, n, nn);
159
160 ri_multiply(wr, wi, fr, fi, n);
161
162 negate(wi, nn);
163 fft_complex_convolution(wr, wi, fr, fi, ldnn);
164
165 make_fft_fract_chirp(wr, wi, v, n, nn);
166
167 ri_multiply(fr, fi, wr, wi, nn);
168
169
170 acopy(wr+n, x, n);
171 acopy(wi+n, y, n);
172
173 delete [] wr;
174 delete [] fr;
175 }
176 // -------------------------
177