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
|
"""
Simpson's rule integration on function `f` with bounds `a` to `b` with `n` intervals.
"""
function integrate(f, a, b, n)
h = (b-a)/2n
x(k) = a + k*h
y(k) = f(x(k))
s = y(0) + y(2n)
for i in 1:2n-1
i % 2 == 1 ? s += 4y(i) : s += 2y(i)
end
I = h/3*s
return I
end
"""
Calculate Fourier coefficient for term with frequency `n` of function `f`.
"""
function c(f, n)
g(t)=f(t)*exp(-2π*im*n*t)
integrate(g,0,1,500)
end
"""
Calculate `n` coefficients for the terms in the fourier series approximation of a given function `f(t)` over domain [0,1].
Corresponding frequencies for the coefficients are 0, 1, -1, 2, -2, and so on.
"""
function calculateCoefficients(f, n)
coefficients = zeros(Complex{Float64}, n)
coefficients[1] = c(f, 0)
Threads.@threads for i in 2:n
coefficients[i] = c(f, (i÷2)*(-1)^i)
end
return coefficients
end
"""
Calculate Fourier series at time `t` with `coefficients` as given by `calculateCoefficients`.
"""
function fourierSeries(t, coefficients)
z = coefficients[1]
n = length(coefficients)
for i in 2:n
z += coefficients[i]*exp(2π*im*(i÷2)*(-1)^i*t)
end
return z
end
"""
Similar to `fourierSeries` but returns a vector of all the intermediate sums. For use in animations.
"""
function fourierArm(t, coefficients)
n = length(coefficients)
arm = zeros(Complex{Float64}, n)
z = 0
for i in 1:n
z += coefficients[i] * exp(2π*im*(i÷2)*(-1)^i*t)
arm[i] = z
end
return arm
end
|
