RosettaCodeData/Task/Fast-Fourier-transform/C-sharp/fast-fourier-transform.cs
2023-07-01 13:44:08 -04:00

103 lines
3.5 KiB
C#

using System;
using System.Numerics;
using System.Linq;
using System.Diagnostics;
// Fast Fourier Transform in C#
public class Program {
/* Performs a Bit Reversal Algorithm on a postive integer
* for given number of bits
* e.g. 011 with 3 bits is reversed to 110 */
public static int BitReverse(int n, int bits) {
int reversedN = n;
int count = bits - 1;
n >>= 1;
while (n > 0) {
reversedN = (reversedN << 1) | (n & 1);
count--;
n >>= 1;
}
return ((reversedN << count) & ((1 << bits) - 1));
}
/* Uses Cooley-Tukey iterative in-place algorithm with radix-2 DIT case
* assumes no of points provided are a power of 2 */
public static void FFT(Complex[] buffer) {
#if false
int bits = (int)Math.Log(buffer.Length, 2);
for (int j = 1; j < buffer.Length / 2; j++) {
int swapPos = BitReverse(j, bits);
var temp = buffer[j];
buffer[j] = buffer[swapPos];
buffer[swapPos] = temp;
}
// Said Zandian
// The above section of the code is incorrect and does not work correctly and has two bugs.
// BUG 1
// The bug is that when you reach and index that was swapped previously it does swap it again
// Ex. binary value n = 0010 and Bits = 4 as input to BitReverse routine and returns 4. The code section above // swaps it. Cells 2 and 4 are swapped. just fine.
// now binary value n = 0010 and Bits = 4 as input to BitReverse routine and returns 2. The code Section
// swap it. Cells 4 and 2 are swapped. WROOOOONG
//
// Bug 2
// The code works on the half section of the cells. In the case of Bits = 4 it means that we are having 16 cells
// The code works on half the cells for (int j = 1; j < buffer.Length / 2; j++) buffer.Length returns 16
// and divide by 2 makes 8, so j goes from 1 to 7. This covers almost everything but what happened to 1011 value
// which must be swap with 1101. and this is the second bug.
//
// use the following corrected section of the code. I have seen this bug in other languages that uses bit
// reversal routine.
#else
for (int j = 1; j < buffer.Length; j++)
{
int swapPos = BitReverse(j, bits);
if (swapPos <= j)
{
continue;
}
var temp = buffer[j];
buffer[j] = buffer[swapPos];
buffer[swapPos] = temp;
}
// First the full length is used and 1011 value is swapped with 1101. Second if new swapPos is less than j
// then it means that swap was happen when j was the swapPos.
#endif
for (int N = 2; N <= buffer.Length; N <<= 1) {
for (int i = 0; i < buffer.Length; i += N) {
for (int k = 0; k < N / 2; k++) {
int evenIndex = i + k;
int oddIndex = i + k + (N / 2);
var even = buffer[evenIndex];
var odd = buffer[oddIndex];
double term = -2 * Math.PI * k / (double)N;
Complex exp = new Complex(Math.Cos(term), Math.Sin(term)) * odd;
buffer[evenIndex] = even + exp;
buffer[oddIndex] = even - exp;
}
}
}
}
public static void Main(string[] args) {
Complex[] input = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0};
FFT(input);
Console.WriteLine("Results:");
foreach (Complex c in input) {
Console.WriteLine(c);
}
}
}