Skip to content

Commit ba2edb3

Browse files
swhardenswharden
authored
Merge pull request #80 from swharden/77
Experimental: Bluesteins algorithm
2 parents 76cd213 + 61fe509 commit ba2edb3

File tree

3 files changed

+194
-0
lines changed

3 files changed

+194
-0
lines changed

‎src/FftSharp.Tests/ExperimentalTests.cs

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,4 +25,24 @@ public void Test_Experimental_DFT()
2525
Assert.AreEqual(expectedImag[i], result[i].Imaginary, 1e-5);
2626
}
2727
}
28+
29+
[Test]
30+
[System.Obsolete]
31+
public void Test_Experimental_Bluestein()
32+
{
33+
double[] prime = { 1, 2, 3, 4, 5, 6, 7 };
34+
System.Numerics.Complex[] result = FftSharp.Experimental.Bluestein(prime);
35+
36+
// tested with python
37+
double[] expectedReal = { 28.00000, -3.50000, -3.50000, -3.50000, -3.50000, -3.50000, -3.50000 };
38+
double[] expectedImag = { -0.00000, 7.26782, 2.79116, 0.79885, -0.79885, -2.79116, -7.26782 };
39+
40+
Assert.AreEqual(expectedReal.Length, result.Length);
41+
42+
for (int i = 0; i < result.Length; i++)
43+
{
44+
Assert.AreEqual(expectedReal[i], result[i].Real, 1e-5);
45+
Assert.AreEqual(expectedImag[i], result[i].Imaginary, 1e-5);
46+
}
47+
}
2848
}

‎src/FftSharp/BluesteinOperations.cs

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
using System;
2+
3+
namespace FftSharp;
4+
5+
internal static class BluesteinOperations
6+
{
7+
/*
8+
* Computes the discrete Fourier transform (DFT) or inverse transform of the given complex vector, storing the result back into the vector.
9+
* The vector can have any length. This is a wrapper function. The inverse transform does not perform scaling, so it is not a true inverse.
10+
*/
11+
public static void Transform(System.Numerics.Complex[] vec, bool inverse)
12+
{
13+
int n = vec.Length;
14+
if (n == 0)
15+
return;
16+
else if ((n & (n - 1)) == 0) // Is power of 2
17+
TransformRadix2(vec, inverse);
18+
else // More complicated algorithm for arbitrary sizes
19+
TransformBluestein(vec, inverse);
20+
}
21+
22+
23+
/*
24+
* Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
25+
* The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
26+
*/
27+
public static void TransformRadix2(System.Numerics.Complex[] vec, bool inverse)
28+
{
29+
// Length variables
30+
int n = vec.Length;
31+
int levels = 0; // compute levels = floor(log2(n))
32+
for (int temp = n; temp > 1; temp >>= 1)
33+
levels++;
34+
if (1 << levels != n)
35+
throw new ArgumentException("Length is not a power of 2");
36+
37+
// Trigonometric table
38+
System.Numerics.Complex[] expTable = new System.Numerics.Complex[n / 2];
39+
double coef = 2 * Math.PI / n * (inverse ? 1 : -1);
40+
for (int i = 0; i < n / 2; i++)
41+
expTable[i] = System.Numerics.Complex.FromPolarCoordinates(1, i * coef);
42+
43+
// Bit-reversed addressing permutation
44+
for (int i = 0; i < n; i++)
45+
{
46+
int j = ReverseBits(i, levels);
47+
if (j > i)
48+
{
49+
System.Numerics.Complex temp = vec[i];
50+
vec[i] = vec[j];
51+
vec[j] = temp;
52+
}
53+
}
54+
55+
// Cooley-Tukey decimation-in-time radix-2 FFT
56+
for (int size = 2; size <= n; size *= 2)
57+
{
58+
int halfsize = size / 2;
59+
int tablestep = n / size;
60+
for (int i = 0; i < n; i += size)
61+
{
62+
for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep)
63+
{
64+
System.Numerics.Complex temp = vec[j + halfsize] * expTable[k];
65+
vec[j + halfsize] = vec[j] - temp;
66+
vec[j] += temp;
67+
}
68+
}
69+
if (size == n) // Prevent overflow in 'size *= 2'
70+
break;
71+
}
72+
}
73+
74+
75+
/*
76+
* Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
77+
* The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
78+
* Uses Bluestein's chirp z-transform algorithm.
79+
*/
80+
public static void TransformBluestein(System.Numerics.Complex[] vec, bool inverse)
81+
{
82+
// Find a power-of-2 convolution length m such that m >= n * 2 + 1
83+
int n = vec.Length;
84+
if (n >= 0x20000000)
85+
throw new ArgumentException("Array too large");
86+
int m = 1;
87+
while (m < n * 2 + 1)
88+
m *= 2;
89+
90+
// Trigonometric table
91+
System.Numerics.Complex[] expTable = new System.Numerics.Complex[n];
92+
double coef = Math.PI / n * (inverse ? 1 : -1);
93+
for (int i = 0; i < n; i++)
94+
{
95+
int j = (int)((long)i * i % (n * 2)); // This is more accurate than j = i * i
96+
expTable[i] = System.Numerics.Complex.Exp(new System.Numerics.Complex(0, j * coef));
97+
}
98+
99+
// Temporary vectors and preprocessing
100+
System.Numerics.Complex[] avec = new System.Numerics.Complex[m];
101+
for (int i = 0; i < n; i++)
102+
avec[i] = vec[i] * expTable[i];
103+
System.Numerics.Complex[] bvec = new System.Numerics.Complex[m];
104+
bvec[0] = expTable[0];
105+
for (int i = 1; i < n; i++)
106+
bvec[i] = bvec[m - i] = System.Numerics.Complex.Conjugate(expTable[i]);
107+
108+
// Convolution
109+
System.Numerics.Complex[] cvec = new System.Numerics.Complex[m];
110+
Convolve(avec, bvec, cvec);
111+
112+
// Postprocessing
113+
for (int i = 0; i < n; i++)
114+
vec[i] = cvec[i] * expTable[i];
115+
}
116+
117+
118+
/*
119+
* Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
120+
*/
121+
public static void Convolve(System.Numerics.Complex[] xvec, System.Numerics.Complex[] yvec, System.Numerics.Complex[] outvec)
122+
{
123+
int n = xvec.Length;
124+
if (n != yvec.Length || n != outvec.Length)
125+
throw new ArgumentException("Mismatched lengths");
126+
xvec = (System.Numerics.Complex[])xvec.Clone();
127+
yvec = (System.Numerics.Complex[])yvec.Clone();
128+
Transform(xvec, false);
129+
Transform(yvec, false);
130+
for (int i = 0; i < n; i++)
131+
xvec[i] *= yvec[i];
132+
Transform(xvec, true);
133+
for (int i = 0; i < n; i++) // Scaling (because this FFT implementation omits it)
134+
outvec[i] = xvec[i] / n;
135+
}
136+
137+
138+
private static int ReverseBits(int val, int width)
139+
{
140+
int result = 0;
141+
for (int i = 0; i < width; i++, val >>= 1)
142+
result = (result << 1) | (val & 1);
143+
return result;
144+
}
145+
}

‎src/FftSharp/Experimental.cs

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,5 +70,34 @@ public static System.Numerics.Complex[] DFT(double[] input, bool inverse = false
7070

7171
return DFT(inputComplex, inverse);
7272
}
73+
74+
/// <summary>
75+
/// Computes the discrete Fourier transform (DFT) of the given array.
76+
/// The array may have any length.
77+
/// Uses Bluestein's chirp z-transform algorithm.
78+
/// </summary>
79+
public static System.Numerics.Complex[] Bluestein(double[] real, bool inverse = false)
80+
{
81+
System.Numerics.Complex[] buffer = new System.Numerics.Complex[real.Length];
82+
83+
for (int i = 0; i < real.Length; i++)
84+
{
85+
buffer[i] = new System.Numerics.Complex(real[i], 0);
86+
}
87+
88+
Bluestein(buffer, inverse);
89+
90+
return buffer;
91+
}
92+
93+
/// <summary>
94+
/// Computes the discrete Fourier transform (DFT) of the given complex vector in place.
95+
/// The vector can have any length.
96+
/// Uses Bluestein's chirp z-transform algorithm.
97+
/// </summary>
98+
public static void Bluestein(System.Numerics.Complex[] buffer, bool inverse = false)
99+
{
100+
BluesteinOperations.TransformBluestein(buffer, inverse);
101+
}
73102
}
74103
}

0 commit comments

Comments
 (0)