SHARE
TWEET

Untitled

a guest Oct 23rd, 2019 99 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. using Plots;
  2. using QuadGK;
  3. using DelimitedFiles;
  4. using Interpolations;
  5.  
  6. struct FourierSeries
  7.     f::Function
  8.     a0::Float64
  9.     ans::Array{Float64,1}
  10.     bns::Array{Float64,1}
  11. end
  12.  
  13. """
  14.     FourierSeries(f; N=10)
  15.  
  16. Generate a FourierSeries object containing the the first `N` (excluding a0)
  17. Fourier coefficients of the function `f` over the interval (-π, π).
  18. """
  19. function FourierSeries(f::Function; N::Int = 10)
  20.     # Compute ans
  21.     a0::Float64 = quadgk(t -> f(t), -π, π)[1] / π
  22.     ans = Array{Float64,1}(undef, N)
  23.     bns = Array{Float64,1}(undef, N)
  24.     for n = 1:N
  25.         ans[n] = quadgk(t -> f(t) * cos(n * t), -π, π)[1] / π
  26.         bns[n] = quadgk(t -> f(t) * sin(n * t), -π, π)[1] / π
  27.     end
  28.     FourierSeries(f, a0, ans, bns)
  29. end
  30.  
  31. """
  32.   eval(t::Float64, fs::FourierSeries)
  33.  
  34. Evaluate `FourierSeries` object at time `t`.
  35. """
  36. function eval(t::Float64, fs::FourierSeries)
  37.     fs.a0 / 2 +
  38.     sum(fs.ans[n] * cos(n * t) + fs.bns[n] * sin(n * t) for n = 1:length(fs.ans))
  39. end
  40.  
  41. """
  42.   eval(ts::AbstractRange, fs::FourierSeries)
  43.  
  44. Evaluate `FourierSeries` object at times `ts`.
  45. """
  46. function eval(ts::AbstractRange, fs::FourierSeries)
  47.     [eval(t, fs) for t in ts]
  48. end
  49.  
  50. """
  51.   circle(x,y,r[; npts=500])
  52.  
  53. Compute x and y coordinates of a circle of radius `r` centered at (`x`, `y`).
  54. """
  55. function circle(x, y, r; npts=500)
  56.     θ = LinRange(0.0, 2π, npts)
  57.     x .+ r * sin.(θ), y .+ r * cos.(θ)
  58. end
  59.  
  60. """
  61.   saw_tooth(t)
  62.  
  63. Returns values of the saw tooth waveform
  64. """
  65. saw_tooth(t) = t
  66.  
  67. """
  68.   square_wave(t)
  69.  
  70. Returns values of the square_wave waveform
  71. """
  72. square_wave(t) = t < 0 ? 0.0 : π
  73.  
  74. """
  75. Make the plot!
  76. """
  77. fs = FourierSeries(saw_tooth)
  78. @gif for i = 1:1500
  79.     shift = i / 100
  80.     plot(
  81.         ts,
  82.         [eval(t + shift, fs) for t in ts],
  83.         ylims = (-1.5, π + 1.5),
  84.         legend = false,
  85.         xlims = (-π, 4π),
  86.         framestyle = :box,
  87.         aspect_ratio = 1
  88.     )
  89.     x = 2.5π
  90.     y = fs.a0 / 2
  91.     r = sqrt(fs.ans[1]^2 + fs.bns[1]^2)
  92.     ϕ = -atan(fs.bns[1], fs.ans[1])
  93.     plot!(
  94.         circle(x, y, r),
  95.         seriestype = [:shape,],
  96.         linecolor = :black,
  97.         c = :white,
  98.         legend = false
  99.     )
  100.     for n = 2:length(fs.ans)
  101.         x += r * sin((n - 1) * (ts[end] + shift) + ϕ)
  102.         y += r * cos((n - 1) * (ts[end] + shift) + ϕ)
  103.         r = sqrt(fs.ans[n]^2 + fs.bns[n]^2)
  104.         ϕ = -atan(fs.bns[n], fs.ans[n])
  105.         plot!(
  106.             circle(x, y, r),
  107.             seriestype = [:shape,],
  108.             linecolor = :black,
  109.             c = :white,
  110.             legend = false,
  111.             fillalpha = 0.0
  112.         )
  113.         scatter!([x], [y], markersize = 0.3)
  114.     end
  115.     x += r * sin(length(fs.ans) * (ts[end] + shift) + ϕ)
  116.     y += r * cos(length(fs.ans) * (ts[end] + shift) + ϕ)
  117.     scatter!([x], [y], markersize = 0.1)
  118.     interps = range(ts[end], stop = x, length = 50)
  119.     plot!(interps, [y for _ in interps], color = :purple)
  120. end every 10
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top