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
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct
{
double re;
double im;
} complex_t;
static complex_t
complex_mul (complex_t a, complex_t b)
{
complex_t r;
r.re = a.re * b.re - a.im * b.im;
r.im = a.re * b.im + a.im * b.re;
return r;
}
static complex_t
complex_add (complex_t a, complex_t b)
{
complex_t r;
r.re = a.re + b.re;
r.im = a.im + b.im;
return r;
}
static complex_t
complex_sub (complex_t a, complex_t b)
{
complex_t r;
r.re = a.re - b.re;
r.im = a.im - b.im;
return r;
}
static void
fft (complex_t *buffer, int n)
{
complex_t *temp, *even, *odd;
complex_t principal, omega;
int i;
if (n == 1)
return;
temp = malloc (n * sizeof (complex_t));
even = temp;
odd = temp + n/2;
for (i = 0; i < n/2; ++i)
{
even[i] = buffer[2 * i];
odd[i] = buffer[2 * i + 1];
}
fft (even, n/2);
fft (odd, n/2);
principal.re = cos (2 * M_PI / n);
principal.im = sin (2 * M_PI / n);
omega.re = 1.0;
omega.im = 0;
for (i = 0; i < n/2; ++i)
{
buffer[i] = complex_add (even[i], complex_mul (omega, odd[i]));
buffer[i + n/2] = complex_sub (even[i], complex_mul (omega, odd[i]));
omega = complex_mul (omega, principal);
}
free (temp);
}
#define N_SAMPLES 128
int
main ()
{
complex_t test[N_SAMPLES];
int i;
for (i = 0; i < N_SAMPLES; ++i)
{
double t = -10 + i * (20 / (double)N_SAMPLES);
double v = (t == 0.0) ? 1.0 : sin (5 * t * M_PI) / (5 * t * M_PI);
test[i].re = v;
test[i].im = 0.0;
}
fft (test, N_SAMPLES);
for (i = 0; i < N_SAMPLES; ++i)
{
printf ("%f %f\n", test[i].re, test[i].im);
}
return 0;
}
|
