Added main scripts and example graphic

1 files changed, 70 insertions, 0 deletions

diff --git a/fourier.jl b/fourier.jl
new file mode 100644
index 0000000..911c984
--- /dev/null
+++ b/fourier.jl
@@ -0,0 +1,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` with 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
\ No newline at end of file