Skip to content

Commit 388799c

Browse files
committed
Merge remote-tracking branch 'upstream/master'
2 parents eb3c7b4 + 4f270e8 commit 388799c

File tree

22 files changed

+635
-367
lines changed

22 files changed

+635
-367
lines changed

README.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
# The Arcane Algorithm Archive
2-
The Arcane Algorithm Archive for all our adventures on LeiosOS / simuleios <br/>
2+
The Arcane Algorithm Archive is a collaborative effort to create a guide for all important algorithms in all languages.
3+
This goal is obviously too ambitious for a book of any size, but it is a great project to learn from and work on and will hopefully become an incredible resource for programmers in the future.
34
The book can be found here: https://www.gitbook.com/book/leios/algorithm-archive/details.
45
The github repository can be found here: https://github.com/leios/algorithm-archive.
56
Most algorithms have been covered on the youtube channel LeiosOS: https://www.youtube.com/user/LeiosOS

book.json

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,12 @@
7171
{
7272
"lang": "elm",
7373
"name": "Elm"
74+
},
75+
{
76+
"lang": "LabVIEW",
77+
"name": "LabVIEW"
7478
}
79+
7580
],
7681
"split": false
7782
}

chapters/FFT/code/c/fft.c

Lines changed: 54 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -1,105 +1,88 @@
1-
// written by Gathros.
2-
31
#include <complex.h>
42
#include <math.h>
5-
6-
// These headers are for presentation not for the algorithm.
73
#include <stdlib.h>
84
#include <time.h>
95
#include <stdio.h>
106

117
#define PI 3.1415926535897932384626
128

13-
void cooley_tukey(double complex *X, const size_t N){
14-
if(N >= 2){
15-
// Splits the array, so the top half are the odd elements and the bottom half are the even ones.
16-
double complex tmp [N/2];
17-
for(size_t i = 0; i < N/2; ++i){
18-
tmp[i] = X[2*i + 1];
19-
X[i] = X[2*i];
20-
}
21-
for(size_t i = 0; i < N/2; ++i){
22-
X[i + N/2] = tmp[i];
23-
}
9+
void cooley_tukey(double complex *X, const size_t N) {
10+
if (N >= 2) {
11+
double complex tmp [N / 2];
12+
for (size_t i = 0; i < N / 2; ++i) {
13+
tmp[i] = X[2*i + 1];
14+
X[i] = X[2*i];
15+
}
16+
for (size_t i = 0; i < N / 2; ++i) {
17+
X[i + N / 2] = tmp[i];
18+
}
2419

25-
// Recursion.
26-
cooley_tukey(X, N/2);
27-
cooley_tukey(X + N/2, N/2);
20+
cooley_tukey(X, N / 2);
21+
cooley_tukey(X + N / 2, N / 2);
2822

29-
// Combine.
30-
for(size_t i = 0; i < N/2; ++i){
31-
X[i + N/2] = X[i] - cexp(-2.0*I*PI*i/N)*X[i + N/2];
32-
X[i] -= (X[i + N/2]-X[i]);
33-
}
34-
}
23+
for (size_t i = 0; i < N / 2; ++i) {
24+
X[i + N / 2] = X[i] - cexp(-2.0 * I * PI * i / N) * X[i + N / 2];
25+
X[i] -= (X[i + N / 2]-X[i]);
26+
}
27+
}
3528
}
3629

37-
void bit_reverse(double complex *X, size_t N){
38-
// Bit reverses the array X[] but only if the size of the array is less then 2^32.
39-
double complex temp;
30+
void bit_reverse(double complex *X, size_t N) {
31+
double complex temp;
4032
unsigned int b;
4133

42-
for(unsigned int i = 0; i < N; ++i){
34+
for (unsigned int i = 0; i < N; ++i) {
4335
b = i;
4436
b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
4537
b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
4638
b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
4739
b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
48-
b = ((b >> 16) | (b << 16)) >> (32 - (unsigned int) log2((double)N));
49-
if(b > i){
40+
b = ((b >> 16) | (b << 16)) >>
41+
(32 - (unsigned int) log2((double)N));
42+
if (b > i) {
5043
temp = X[b];
5144
X[b] = X[i];
5245
X[i] = temp;
5346
}
5447
}
5548
}
5649

57-
void iterative_cooley_tukey(double complex *X, size_t N){
58-
int stride;
59-
double complex v,w;
60-
61-
// Bit reverse the array.
62-
bit_reverse(X, N);
63-
64-
// Preform the butterfly on the array.
65-
for(int i = 1; i <= log2((double)N); ++i){
66-
stride = pow(2, i);
67-
w = cexp(-2.0*I*PI/stride);
68-
for(size_t j = 0; j < N; j += stride){
69-
v = 1.0;
70-
for(size_t k = 0; k < stride/2; k++){
71-
X[k + j + stride/2] = X[k + j] - v*X[k + j + stride/2];
72-
X[k + j] -= (X[k + j + stride/2] - X[k + j]);
73-
v *= w;
74-
}
75-
}
76-
}
50+
void iterative_cooley_tukey(double complex *X, size_t N) {
51+
bit_reverse(X, N);
52+
53+
for (int i = 1; i <= log2((double)N); ++i) {
54+
int stride = pow(2, i);
55+
double complex w = cexp(-2.0 * I * PI / stride);
56+
for (size_t j = 0; j < N; j += stride) {
57+
double complex v = 1.0;
58+
for (size_t k = 0; k < stride / 2; ++k) {
59+
X[k + j + stride / 2] = X[k + j] - v * X[k + j + stride / 2];
60+
X[k + j] -= (X[k + j + stride / 2] - X[k + j]);
61+
v *= w;
62+
}
63+
}
64+
}
7765
}
7866

79-
void approx(double complex *X, double complex *Y, size_t N){
80-
// This is to show that the arrays are approximate.
81-
for(size_t i = 0; i < N; ++i){
82-
printf("%f\n", cabs(X[i]) - cabs(Y[i]));
83-
}
67+
void approx(double complex *X, double complex *Y, size_t N) {
68+
for (size_t i = 0; i < N; ++i) {
69+
printf("%f\n", cabs(X[i]) - cabs(Y[i]));
70+
}
8471
}
8572

86-
int main(){
87-
// Initalizing the arrays for FFT.
88-
srand(time(NULL));
89-
const size_t N = 64;
90-
double complex x[N], y[N], z[N];
91-
for(size_t i = 0; i < N; ++i){
92-
x[i] = rand() / (double) RAND_MAX;
93-
y[i] = x[i];
94-
z[i] = x[i];
95-
}
73+
int main() {
74+
srand(time(NULL));
75+
double complex x[64], y[64], z[64];
76+
for (size_t i = 0; i < 64; ++i) {
77+
x[i] = rand() / (double) RAND_MAX;
78+
y[i] = x[i];
79+
z[i] = x[i];
80+
}
9681

97-
// Preform FFT.
98-
cooley_tukey(y, N);
99-
iterative_cooley_tukey(z, N);
82+
cooley_tukey(y, 64);
83+
iterative_cooley_tukey(z, 64);
10084

101-
// Check if the different methods are approximate.
102-
approx(y, z, N);
85+
approx(y, z, 64);
10386

104-
return 0;
87+
return 0;
10588
}

chapters/FFT/cooley_tukey.md

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,9 @@ If we take a sum sinusoidal functions (like $$\sin(\omega t)$$ or $$\cos(\omega
5151
Each constituent wave can be described by only one value: $$\omega$$.
5252
So, instead of representing these curves as seen above, we could instead describe them as peaks in frequency space, as shown below.
5353

54-
![Fourier Example](res/FT_example.png)
54+
<p align="center">
55+
<img src="res/FT_example.png" width="500" height="250" />
56+
</p>
5557

5658
This is what the Fourier Transform does!
5759
After performing the transform, it is now much, much easier to understand precisely which frequencies are in our waveform, which is essential to most areas of signal processing.
@@ -138,9 +140,9 @@ In the end, the code looks like:
138140
{% sample lang="jl" %}
139141
[import:14-31, lang:"julia"](code/julia/fft.jl)
140142
{% sample lang="c" %}
141-
[import:13-35, lang:"c_cpp"](code/c/fft.c)
143+
[import:9-28, lang:"c_cpp"](code/c/fft.c)
142144
{% sample lang="cpp" %}
143-
[import:19-44, lang:"c_cpp"](code/c++/fft.cpp)
145+
[import:27-57, lang:"c_cpp"](code/c++/fft.cpp)
144146
{% sample lang="hs" %}
145147
[import:6-19, lang:"haskell"](code/hs/fft.hs)
146148
{% sample lang="py2" %}
@@ -170,7 +172,9 @@ And at each step, we use the appropriate term.
170172
For example, imagine we need to perform an FFT of an array of only 2 elements.
171173
We can represent this addition with the following (radix-2) butterfly:
172174

173-
![Radix-2, positive W](res/radix-2screen_positive.jpg)
175+
<p align="center">
176+
<img src="res/radix-2screen_positive.jpg" width="400" height="225" />
177+
</p>
174178

175179
Here, the diagram means the following:
176180

@@ -182,7 +186,9 @@ $$
182186

183187
However, it turns out that the second half of our array of $$\omega$$ values is always the negative of the first half, so $$\omega_2^0 = -\omega_2^1$$, so we can use the following butterfly diagram:
184188

185-
![Radix-2](res/radix-2screen.jpg)
189+
<p align="center">
190+
<img src="res/radix-2screen.jpg" width="400" />
191+
</p>
186192

187193
With the following equations:
188194

@@ -197,14 +203,18 @@ Now imagine we need to combine more elements.
197203
In this case, we start with simple butterflies, as shown above, and then sum butterflies of butterflies.
198204
For example, if we have 8 elements, this might look like this:
199205

200-
![Radix-8](res/radix-8screen.jpg)
206+
<p align="center">
207+
<img src="res/radix-8screen.jpg" width="500" height="500" />
208+
</p>
201209

202210
Note that we can perform a DFT directly before using any butterflies, if we so desire, but we need to be careful with how we shuffle our array if that's the case.
203211
In the code snippet provided in the previous section, the subdivision was performed in the same function as the concatenation, so the ordering was always correct; however, if we were to re-order with bit-reversal, this might not be the case.
204212

205213
For example, take a look at the ordering of FFT ([found on wikipedia](https://en.wikipedia.org/wiki/Butterfly_diagram)) that performs the DFT shortcut:
206214

207-
![Butterfly Diagram](res/butterfly_diagram.png)
215+
<p align="center">
216+
<img src="res/butterfly_diagram.png" width="600" height="500" />
217+
</p>
208218

209219
Here, the ordering of the array was simply divided into even and odd elements once, but they did not recursively divide the arrays of even and odd elements again because they knew they would perform a DFT soon thereafter.
210220

chapters/convolutions/code/c/convolutions.c

Lines changed: 35 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -5,50 +5,49 @@
55

66
#include <math.h>
77

8-
void fft(double complex *X, size_t N){
9-
if(N >= 2){
10-
// Splits the array, so the top half are the odd elements and the bottom half are the even ones.
11-
double complex tmp [N/2];
12-
for(size_t i = 0; i < N/2; ++i){
13-
tmp[i] = X[2*i + 1];
14-
X[i] = X[2*i];
8+
void fft(double complex *X, size_t N) {
9+
if (N >= 2) {
10+
double complex tmp [N / 2];
11+
for (size_t i = 0; i < N / 2; ++i) {
12+
tmp[i] = X[2 * i + 1];
13+
X[i] = X[2 * i];
1514
}
16-
for(size_t i = 0; i < N/2; ++i){
17-
X[i + N/2] = tmp[i];
15+
16+
for (size_t i = 0; i < N / 2; ++i) {
17+
X[i + N / 2] = tmp[i];
1818
}
1919

20-
// Recursion.
21-
fft(X, N/2);
22-
fft(X + N/2, N/2);
20+
fft(X, N / 2);
21+
fft(X + N / 2, N / 2);
2322

24-
// Combine.
25-
for(size_t i = 0; i < N/2; ++i){
26-
X[i + N/2] = X[i] - cexp(-2.0*I*M_PI*i/N)*X[i + N/2];
27-
X[i] -= (X[i + N/2]-X[i]);
23+
for (size_t i = 0; i < N / 2; ++i) {
24+
X[i + N/2] = X[i] - cexp(-2.0 * I * M_PI * i / N) * X[i + N / 2];
25+
X[i] -= (X[i + N / 2] - X[i]);
2826
}
2927
}
3028
}
3129

32-
void ifft(double complex *x, size_t n){
33-
for(size_t i = 0; i < n; i++){
30+
void ifft(double complex *x, size_t n) {
31+
for (size_t i = 0; i < n; ++i) {
3432
x[i] = conj(x[i]);
3533
}
3634

3735
fft(x, n);
3836

39-
for(size_t i = 0; i < n; i++){
40-
x[i] = conj(x[i])/n;
37+
for (size_t i = 0; i < n; ++i) {
38+
x[i] = conj(x[i]) / n;
4139
}
4240
}
4341

4442
// This section is a part of the algorithm
4543

46-
void conv(double complex *signal1, double complex *signal2, double complex* out, size_t n1, size_t n2){
44+
void conv(double complex *signal1, double complex *signal2, double complex* out,
45+
size_t n1, size_t n2) {
4746
double complex sum = 0;
4847

49-
for(size_t i = 0; i < (n1 < n2? n2 : n1); i++){
50-
for(size_t j = 0; j < i; j++){
51-
if(j < n1){
48+
for (size_t i = 0; i < (n1 < n2? n2 : n1); ++i) {
49+
for (size_t j = 0; j < i; ++j) {
50+
if (j < n1) {
5251
sum += signal1[j] * signal2[i-j];
5352
}
5453
}
@@ -57,22 +56,24 @@ void conv(double complex *signal1, double complex *signal2, double complex* out,
5756
}
5857
}
5958

60-
void conv_fft(double complex *signal1, double complex *signal2, double complex* out, size_t n){
59+
void conv_fft(double complex *signal1, double complex *signal2,
60+
double complex* out, size_t n) {
6161
fft(signal1, n);
6262
fft(signal2, n);
6363

64-
for(size_t i = 0; i < n; i++){
65-
out[i] = signal1[i]*signal2[i];
64+
for (size_t i = 0; i < n; ++i) {
65+
out[i] = signal1[i] * signal2[i];
6666
}
6767

6868
ifft(out, n);
6969
}
7070

71-
int main(){
72-
double complex signal1[64], signal2[64], signal3[128], signal4[128], out1[128], out2[128];
71+
int main() {
72+
double complex signal1[64], signal2[64], signal3[128], signal4[128],
73+
out1[128], out2[128];
7374

74-
for(size_t i = 0; i < 128; i++){
75-
if(i >= 16 && i < 48){
75+
for (size_t i = 0; i < 128; ++i) {
76+
if (i >= 16 && i < 48) {
7677
signal1[i] = 1.0;
7778
signal2[i] = 1.0;
7879
signal3[i] = 1.0;
@@ -97,8 +98,9 @@ int main(){
9798
conv(signal1, signal2, out1, 64, 64);
9899
conv_fft(signal3, signal4, out2, 128);
99100

100-
for(size_t i = 0; i < 64; i++){
101-
printf("%zu %f %+fi\n", i, creal(out1[i]) - creal(out2[i]), cimag(out1[i]) - cimag(out2[i]));
101+
for (size_t i = 0; i < 64; ++i) {
102+
printf("%zu %f %+fi\n", i, creal(out1[i]) - creal(out2[i]),
103+
cimag(out1[i]) - cimag(out2[i]));
102104
}
103105

104106
return 0;

0 commit comments

Comments
 (0)