repl.it
@wolf32669/

C FFT-3

C

No description

fork
loading
Files
  • main.c
main.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
math.random():
#include <stdlib.h>
#include <math.h>

#define PI 3.14159265f

typedef float num; // change to double for precision

typedef struct complex {
  num r;
  num c;
} complex;

complex c_root(num phase) {
  return (complex){cos(phase), sin(phase)};
}

// 00000004
// 00040000
// 04000000
// 40000000
// 10000000
// 20000000
unsigned int bit_reverse(unsigned int n, int size) {
  n = ((n & 0x0000ffff) << 16) | ((n & 0xffff0000) >> 16);
  n = ((n & 0x00ff00ff) << 8)  | ((n & 0xff00ff00) >> 8);
  n = ((n & 0x0f0f0f0f) << 4)  | ((n & 0xf0f0f0f0) >> 4);
  n = ((n & 0x33333333) << 2)  | ((n & 0xCCCCCCCC) >> 2);
  n = ((n & 0x55555555) << 1)  | ((n & 0xAAAAAAAA) >> 1);
  return n >> (32 - size);
}

void bit_reverse_perm(complex* nums, int size) {
  for (int i = 0; i < (1 << size); i++) {
    int j = bit_reverse(i, size);
    if (i < j) {
      complex temp = nums[i];
      nums[i] = nums[j];
      nums[j] = temp;
    }
  }
}

int get_size(int n) {
  int res = 0;
  while (n > 1) {n <<= 1; res++;}
  return res;
}

// takes argument bit-reversed
void fft(complex* nums, int size) {
  if (size == 0) {
    return;
  }

  complex* left = nums;
  complex* right = nums + (1 << (size - 1));
  fft(left, size - 1);
  fft(right, size - 1);
  num phase = -2 * PI / (1 << size);
  for (int i = 0; i < 1 << (size - 1); i++) {
    complex l = left[i];
    complex r = right[i];
    complex root = c_root(phase * i);
    r = (complex){r.r * root.r - r.c * root.c,
                  r.r * root.c + r.c * root.r};

    left[i] = (complex){l.r + r.r, l.c + r.c};
    right[i] = (complex){l.r - r.r, l.c - r.c};
  }
}

// takes argument bit-reversed
void ifft(complex* nums, int size) {
  int n = 1 << size;
  for (int i = 0; i < n; i++) {
    nums[i] = (complex){nums[i].r, -nums[i].c};
  }
  fft(nums, size);
  for (int i = 0; i < n; i++) {
    nums[i] = (complex){nums[i].r / n, -nums[i].c / n};
  }
}

void print_nums(complex* nums, int len, char* fmt) {
  for (int i = 0; i < len; i++) {
    printf(fmt, nums[i].r, nums[i].c);
  }
}

int main(void) {
  complex* l = malloc(4 * sizeof(complex));

  l[0] = (complex){+1, +0};
  l[1] = (complex){-1, +0};
  l[2] = (complex){+1, +0};
  l[3] = (complex){-1, +0};
  printf("start:\n");
  print_nums(l, 4, "  %f + %fi\n");

  bit_reverse_perm(l, 2);
  fft(l, 2);
  printf("fft:\n");
  print_nums(l, 4, "  %f + %fi\n");
  
  bit_reverse_perm(l, 2);
  ifft(l, 2);
  printf("inverted:\n");
  print_nums(l, 4, "  %f + %fi\n");
  return 0;
}
Fetching token
?