Décomposition de Fourier d'une courbe fermée dans $R^2$

\documentclass{article}

\usepackage{luapackageloader}
\usepackage{luamplib}

\directlua{dofile("Fourier.lua")}


\pagestyle{empty}

\begin{document}


\directlua{
  plotDecompAnim("f",300,50,0.26,360)
}



\end{document}

complex = require "complex"
local mpkpse = kpse.new('luatex', 'mpost')

local function finder(name, mode, ftype)
   if mode == "w" then
  return name
   else
       if ftype == 'pfb' then ftype='type1 fonts' end
       return  mpkpse:find_file(name,ftype)
   end
end


function getpathfrommp(s,nbrPoints,scale)
   local mp = mplib.new({
         find_file = finder,})
   mp:execute('input metafun ;')
   local rettable
   rettable = mp:execute('fontmapfile "=lm-ec.map"; picture lettre; path contourLettre; lettre := glyph "' .. s .. '" of "ec-lmri10"; path p; beginfig(1); for item within lettre: contourLettre := pathpart item; for i:=1 upto'.. nbrPoints ..': if i=1: p := point i/'..nbrPoints..' along contourLettre; else: p:= p--(point i/'..nbrPoints..' along contourLettre);fi; endfor; draw p scaled '..scale..'; endfor; endfig;end;')
   output = {}
   if rettable.status == 0 then
      figures = rettable.fig
      figure = figures[1]
      local objects = figure:objects()
      local segment = objects[1]
      for point =1, #segment.path do
         output[point] = {}
         output[point].x = segment.path[point].x_coord
         output[point].y = segment.path[point].y_coord
      end
   else
      print("error")
   end
   return output
end



function pathToComplex(path)
   local complexPath
   complexPath = {}
   for i=1,#path do
      complexPath[i] = complex.new(path[i].x,path[i].y)
   end
   return complexPath
end

function listComplexToPolar(cplxList)
   local output = {}
   for i=1,#cplxList do
      output[i] = {}
      output[i].abs,output[i].arg = complex.polar(cplxList[i])
   end
   return output
end

function cn(f,n)
   local CN
   CN = complex.new(0.0,0.0)
   local N = #f
   for i=0,N-1 do
      exposant = complex.new(0.0,-2.0*math.pi*n*i/N)
      Exp = complex.exp(exposant)
      CN = complex.add(CN,complex.mulnum(complex.mul(f[i+1],Exp),1.0/N))
   end
   return CN
end

function cnList(f)
   local CNlist = {}
   local N = #f
   for i=0,N do
      --print(math.floor(i-N/2))
      CNlist[i] = cn(f,math.floor(i-N/2))
   end
   return CNlist
end



function drawPath(path)
   local str
   str = "\\begin{mplibcode}\n beginfig(1);\n"
   str = str.."path p;\n"
   str = str.." p:="
   for i=1,#path do
      str = str.."("..string.format("%f",path[i].x)..","..string.format("%f",path[i].y)..")--"
   end
   str = str.."cycle;\n"
   str = str.."draw p scaled 0.1;\n endfig;\n"
   str = str.."\\end{mplibcode}\n"
   --print(str)
   tex.sprint(str)
end

function mpCodePath(path)
   local str
   str = ""
   str = str.."path p;\n"
   str = str.." p:="
   for i=1,#path do
      str = str.."("..string.format("%f",path[i].x)..","..string.format("%f",path[i].y)..")--"
   end
   str = str.."cycle;\n"
   str = str.."draw p;\n"
   str = str.. "pair ll, lr, ur, ul; ll := llcorner p; ur := urcorner p; lr := lrcorner p; ul := ulcorner p;\n"
   str = str.. "Wdth := abs(xpart lr - xpart ll); Hght := abs(ypart ul - ypart ll); prcW := 0.8;prcH := 0.3;\n"
   str = str.. "draw (ll+(-prcW*Wdth,-prcH*Hght))--(lr+(+prcW*Wdth,-prcH*Hght))--(ur+(+prcW*Wdth,+prcH*Hght))--(ul+(-prcW*Wdth,+prcH*Hght))--cycle;\n"
   return str
end

function mpCodeCircle(cn,shift)
   local str
   local abs,arg
   abs,arg= complex.polar(cn)
   str = "draw fullcircle scaled "..string.format("%f",2*abs).." shifted ("..string.format("%f",shift[1])..","..string.format("%f",shift[2])..") withcolor (0.7,0.7,0.7);\n"
   str = str.."drawarrow ((0,0)--("..string.format("%f",cn[1])..","..string.format("%f",cn[2])..") )  shifted ("..string.format("%f",shift[1])..","..string.format("%f",shift[2])..") withpen pencircle scaled 1pt withcolor (0.7,0.3,0.3);"
   return str
end

function coreDecomp(path, complexPath, cnList,nbrFourier,t)
   local str
   local cnListRotated = {}
   local zero = math.floor(#cnList/2)
   local NFourier = math.floor(nbrFourier/2)
   cnListRotated[zero] = cnList[zero]
   for i=1, zero do
      cnListRotated[zero+i] = complex.mul(cnList[zero+i],complex.exp(complex.new(0.0, i*2*math.pi*t)))
      cnListRotated[zero-i] = complex.mul(cnList[zero-i],complex.exp(complex.new(0.0, -i*2*math.pi*t)))
   end
   --local listAbsArg = listComplexToPolar(cnList)
   local str = "\\begin{mplibcode}\nverbatimtex \\leavevmode etex;beginfig(1);"
   local mpCodeLetter = mpCodePath(path)
   str = str..mpCodeLetter
   --for i=1,#cnList do

   local currentC = complex.new(0,0)
   str = str .. mpCodeCircle(cnListRotated[zero],currentC)
   currentC = complex.add(currentC,cnListRotated[zero])
   for i=1,NFourier do
      str = str .. mpCodeCircle(cnListRotated[zero+i],currentC)
      currentC = complex.add(currentC,cnListRotated[zero+i])
      str = str .. mpCodeCircle(cnListRotated[zero-i],currentC)
      currentC = complex.add(currentC,cnListRotated[zero-i])
   end
   --end
   str = str .. "endfig;\n\\end{mplibcode}\n\\newpage"
   --print(str)
   return str
end


function plotDecomp(letter,nbrPoints,nbrFourier,t,scale)
   local str
   local path = getpathfrommp(letter,nbrPoints,scale)
   local complexPath = pathToComplex(path)
   local cnList = cnList(complexPath)
   str = coreDecomp(path,complexPath,cnList,nbrFourier,t)
   tex.sprint(str)
end

function plotDecompAnim(letter,nbrPoints,nbrFourier,scale,nbrFrame)
   local str
   local path = getpathfrommp(letter,nbrPoints,scale)
   local complexPath = pathToComplex(path)
   local cnList = cnList(complexPath)
   local cnListRotated = {}
   local zero = math.floor(#cnList/2)
   local NFourier = math.floor(nbrFourier/2)
   for frame=0,nbrFrame-1 do
      t=frame/nbrFrame
      str = coreDecomp(path,complexPath,cnList,nbrFourier,t)
      tex.sprint(str)
   end
end

-- complex 0.3.0
-- Lua 5.1

-- 'complex' provides common tasks with complex numbers

-- function complex.to( arg ); complex( arg )
-- returns a complex number on success, nil on failure
-- arg := number or { number,number } or ( "(-)" and/or "(+/-)i" )
--      e.g. 5; {2,3}; "2", "2+i", "-2i", "2^2*3+1/3i"
--      note: 'i' is always in the numerator, spaces are not allowed

-- a complex number is defined as carthesic complex number
-- complex number := { real_part, imaginary_part }
-- this gives fast access to both parts of the number for calculation
-- the access is faster than in a hash table
-- the metatable is just a add on, when it comes to speed, one is faster using a direct function call

-- http://luaforge.net/projects/LuaMatrix
-- http://lua-users.org/wiki/ComplexNumbers

-- Licensed under the same terms as Lua itself.

--/////////////--
--// complex //--
--/////////////--

-- link to complex table
local complex = {}

-- link to complex metatable
local complex_meta = {}

-- complex.to( arg )
-- return a complex number on success
-- return nil on failure
local _retone = function() return 1 end
local _retminusone = function() return -1 end
function complex.to( num )
   -- check for table type
   if type( num ) == "table" then
      -- check for a complex number
      if getmetatable( num ) == complex_meta then
         return num
      end
      local real,imag = tonumber( num[1] ),tonumber( num[2] )
      if real and imag then
         return setmetatable( { real,imag }, complex_meta )
      end
      return
   end
   -- check for number
   local isnum = tonumber( num )
   if isnum then
      return setmetatable( { isnum,0 }, complex_meta )
   end
   if type( num ) == "string" then
      -- check for real and complex
      -- number chars [%-%+%*%^%d%./Ee]
      local real,sign,imag = string.match( num, "^([%-%+%*%^%d%./Ee]*%d)([%+%-])([%-%+%*%^%d%./Ee]*)i$" )
      if real then
         if string.lower(string.sub(real,1,1)) == "e"
         or string.lower(string.sub(imag,1,1)) == "e" then
            return
         end
         if imag == "" then
            if sign == "+" then
               imag = _retone
            else
               imag = _retminusone
            end
         elseif sign == "+" then
            imag = loadstring("return tonumber("..imag..")")
         else
            imag = loadstring("return tonumber("..sign..imag..")")
         end
         real = loadstring("return tonumber("..real..")")
         if real and imag then
            return setmetatable( { real(),imag() }, complex_meta )
         end
         return
      end
      -- check for complex
      local imag = string.match( num,"^([%-%+%*%^%d%./Ee]*)i$" )
      if imag then
         if imag == "" then
            return setmetatable( { 0,1 }, complex_meta )
         elseif imag == "-" then
            return setmetatable( { 0,-1 }, complex_meta )
         end
         if string.lower(string.sub(imag,1,1)) ~= "e" then
            imag = loadstring("return tonumber("..imag..")")
            if imag then
               return setmetatable( { 0,imag() }, complex_meta )
            end
         end
         return
      end
      -- should be real
      local real = string.match( num,"^(%-*[%d%.][%-%+%*%^%d%./Ee]*)$" )
      if real then
         real = loadstring( "return tonumber("..real..")" )
         if real then
            return setmetatable( { real(),0 }, complex_meta )
         end
      end
   end
end

-- complex( arg )
-- same as complex.to( arg )
-- set __call behaviour of complex
setmetatable( complex, { __call = function( _,num ) return complex.to( num ) end } )

-- complex.new( real, complex )
-- fast function to get a complex number, not invoking any checks
function complex.new( ... )
   return setmetatable( { ... }, complex_meta )
end

-- complex.type( arg )
-- is argument of type complex
function complex.type( arg )
   if getmetatable( arg ) == complex_meta then
      return "complex"
   end
end

-- complex.convpolar( r, phi )
-- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
-- r (radius) is a number
-- phi (angle) must be in radians; e.g. [0 - 2pi]
function complex.convpolar( radius, phi )
   return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
end

-- complex.convpolardeg( r, phi )
-- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
-- r (radius) is a number
-- phi must be in degrees; e.g. [0° - 360°]
function complex.convpolardeg( radius, phi )
   phi = phi/180 * math.pi
   return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
end

--// complex number functions only

-- complex.tostring( cx [, formatstr] )
-- to string or real number
-- takes a complex number and returns its string value or real number value
function complex.tostring( cx,formatstr )
   local real,imag = cx[1],cx[2]
   if formatstr then
      if imag == 0 then
         return string.format( formatstr, real )
      elseif real == 0 then
         return string.format( formatstr, imag ).."i"
      elseif imag > 0 then
         return string.format( formatstr, real ).."+"..string.format( formatstr, imag ).."i"
      end
      return string.format( formatstr, real )..string.format( formatstr, imag ).."i"
   end
   if imag == 0 then
      return real
   elseif real == 0 then
      return ((imag==1 and "") or (imag==-1 and "-") or imag).."i"
   elseif imag > 0 then
      return real.."+"..(imag==1 and "" or imag).."i"
   end
   return real..(imag==-1 and "-" or imag).."i"
end

-- complex.print( cx [, formatstr] )
-- print a complex number
function complex.print( ... )
   print( complex.tostring( ... ) )
end

-- complex.polar( cx )
-- from complex number to polar coordinates
-- output in radians; [-pi,+pi]
-- returns r (radius), phi (angle)
function complex.polar( cx )
   return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] )
end

-- complex.polardeg( cx )
-- from complex number to polar coordinates
-- output in degrees; [-180°,180°]
-- returns r (radius), phi (angle)
function complex.polardeg( cx )
   return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] ) / math.pi * 180
end

-- complex.mulconjugate( cx )
-- multiply with conjugate, function returning a number
function complex.mulconjugate( cx )
   return cx[1]^2 + cx[2]^2
end

-- complex.abs( cx )
-- get the absolute value of a complex number
function complex.abs( cx )
   return math.sqrt( cx[1]^2 + cx[2]^2 )
end

-- complex.get( cx )
-- returns real_part, imaginary_part
function complex.get( cx )
   return cx[1],cx[2]
end

-- complex.set( cx, real, imag )
-- sets real_part = real and imaginary_part = imag
function complex.set( cx,real,imag )
   cx[1],cx[2] = real,imag
end

-- complex.is( cx, real, imag )
-- returns true if, real_part = real and imaginary_part = imag
-- else returns false
function complex.is( cx,real,imag )
   if cx[1] == real and cx[2] == imag then
      return true
   end
   return false
end

--// functions returning a new complex number

-- complex.copy( cx )
-- copy complex number
function complex.copy( cx )
   return setmetatable( { cx[1],cx[2] }, complex_meta )
end

-- complex.add( cx1, cx2 )
-- add two numbers; cx1 + cx2
function complex.add( cx1,cx2 )
   return setmetatable( { cx1[1]+cx2[1], cx1[2]+cx2[2] }, complex_meta )
end

-- complex.sub( cx1, cx2 )
-- subtract two numbers; cx1 - cx2
function complex.sub( cx1,cx2 )
   return setmetatable( { cx1[1]-cx2[1], cx1[2]-cx2[2] }, complex_meta )
end

-- complex.mul( cx1, cx2 )
-- multiply two numbers; cx1 * cx2
function complex.mul( cx1,cx2 )
   return setmetatable( { cx1[1]*cx2[1] - cx1[2]*cx2[2],cx1[1]*cx2[2] + cx1[2]*cx2[1] }, complex_meta )
end

-- complex.mulnum( cx, num )
-- multiply complex with number; cx1 * num
function complex.mulnum( cx,num )
   return setmetatable( { cx[1]*num,cx[2]*num }, complex_meta )
end

-- complex.div( cx1, cx2 )
-- divide 2 numbers; cx1 / cx2
function complex.div( cx1,cx2 )
   -- get complex value
   local val = cx2[1]^2 + cx2[2]^2
   -- multiply cx1 with conjugate complex of cx2 and divide through val
   return setmetatable( { (cx1[1]*cx2[1]+cx1[2]*cx2[2])/val,(cx1[2]*cx2[1]-cx1[1]*cx2[2])/val }, complex_meta )
end

-- complex.divnum( cx, num )
-- divide through a number
function complex.divnum( cx,num )
   return setmetatable( { cx[1]/num,cx[2]/num }, complex_meta )
end

-- complex.pow( cx, num )
-- get the power of a complex number
function complex.pow( cx,num )
   if math.floor( num ) == num then
      if num < 0 then
         local val = cx[1]^2 + cx[2]^2
         cx = { cx[1]/val,-cx[2]/val }
         num = -num
      end
      local real,imag = cx[1],cx[2]
      for i = 2,num do
         real,imag = real*cx[1] - imag*cx[2],real*cx[2] + imag*cx[1]
      end
      return setmetatable( { real,imag }, complex_meta )
   end
   -- we calculate the polar complex number now
   -- since then we have the versatility to calc any potenz of the complex number
   -- then we convert it back to a carthesic complex number, we loose precision here
   local length,phi = math.sqrt( cx[1]^2 + cx[2]^2 )^num, math.atan2( cx[2], cx[1] )*num
   return setmetatable( { length * math.cos( phi ), length * math.sin( phi ) }, complex_meta )
end

-- complex.sqrt( cx )
-- get the first squareroot of a complex number, more accurate than cx^.5
function complex.sqrt( cx )
   local len = math.sqrt( cx[1]^2+cx[2]^2 )
   local sign = (cx[2]<0 and -1) or 1
   return setmetatable( { math.sqrt((cx[1]+len)/2), sign*math.sqrt((len-cx[1])/2) }, complex_meta )
end

-- complex.ln( cx )
-- natural logarithm of cx
function complex.ln( cx )
   return setmetatable( { math.log(math.sqrt( cx[1]^2 + cx[2]^2 )),
      math.atan2( cx[2], cx[1] ) }, complex_meta )
end

-- complex.exp( cx )
-- exponent of cx (e^cx)
function complex.exp( cx )
   local expreal = math.exp(cx[1])
   return setmetatable( { expreal*math.cos(cx[2]), expreal*math.sin(cx[2]) }, complex_meta )
end

-- complex.conjugate( cx )
-- get conjugate complex of number
function complex.conjugate( cx )
   return setmetatable( { cx[1], -cx[2] }, complex_meta )
end

-- complex.round( cx [,idp] )
-- round complex numbers, by default to 0 decimal points
function complex.round( cx,idp )
   local mult = 10^( idp or 0 )
   return setmetatable( { math.floor( cx[1] * mult + 0.5 ) / mult,
      math.floor( cx[2] * mult + 0.5 ) / mult }, complex_meta )
end

--// metatable functions

complex_meta.__add = function( cx1,cx2 )
   local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
   return complex.add( cx1,cx2 )
end
complex_meta.__sub = function( cx1,cx2 )
   local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
   return complex.sub( cx1,cx2 )
end
complex_meta.__mul = function( cx1,cx2 )
   local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
   return complex.mul( cx1,cx2 )
end
complex_meta.__div = function( cx1,cx2 )
   local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
   return complex.div( cx1,cx2 )
end
complex_meta.__pow = function( cx,num )
   if num == "*" then
      return complex.conjugate( cx )
   end
   return complex.pow( cx,num )
end
complex_meta.__unm = function( cx )
   return setmetatable( { -cx[1], -cx[2] }, complex_meta )
end
complex_meta.__eq = function( cx1,cx2 )
   if cx1[1] == cx2[1] and cx1[2] == cx2[2] then
      return true
   end
   return false
end
complex_meta.__tostring = function( cx )
   return tostring( complex.tostring( cx ) )
end
complex_meta.__concat = function( cx,cx2 )
   return tostring(cx)..tostring(cx2)
end
-- cx( cx, formatstr )
complex_meta.__call = function( ... )
   print( complex.tostring( ... ) )
end
complex_meta.__index = {}
for k,v in pairs( complex ) do
   complex_meta.__index[k] = v
end

return complex

--///////////////--
--// chillcode //--
--///////////////--