your-own-drum/wave.lua

323 lines
7.3 KiB
Lua

--Render and simulate 1D wave equation.
local love = love
local N = 33
local SOUNDSPEED
local IMPULSESIZE
local DAMPING
local old, cur, new --States add beginning of last tick, current tick, and next tick respectively.
local function Current() return cur end
local function Next() return new end
local Integrate
local Interpolate
local Derivative
local SecondDerivative
local DFT
local shader = love.graphics.newShader([[
uniform float re[33];
uniform float im[33];
uniform float score;
//Slow IDFT
float r( float x )
{
float r = re[0];
for( int k = 1; k < 17; k++ )
{
float c = cos( x * float( k ) );
float s = sin( x * float( k ) );
r +=
+ c * re[k]
- s * im[k]
+ c * re[33 - k]
+ s * im[33 - k];
}
return r;
}
vec4 effect(vec4 color, Image tex, vec2 texture_coords, vec2 screen_coords)
{
vec2 p = (3.0 * screen_coords - 1.5 * love_ScreenSize.xy ) / love_ScreenSize.y;
p.y = -p.y;
float r = r( atan(p.y, p.x) ) - length( p );
float q = float( r < 0.05 ) * clamp( 1.0 - score, 0.0, 1.0 ) ;
return vec4( q + (1.0 + clamp( score, 0.0, 1.0 ) * r * r * 0.3) * color.xyz, float(r > 0.0) ) ;
}
]])
--Calculate discrete fourier transform of radius function.
do
local twiddlec = {}
local twiddles = {}
for i = 1, N do
twiddlec[i], twiddles[i] = math.cos( - 2.0 * ( i - 1 ) * math.pi / N ), math.sin( - 2.0 * (i - 1) * math.pi / N )
end
local function Twiddle( j )
j = 1 + j % N
return twiddlec[j], twiddles[j]
end
--Slow Discrete Fourier Transform of a real sequence.
DFT = function( wave )
local seq = wave.radii
local dre, dim = wave.dftre, wave.dftim
for k = 0, N - 1 do
local x, y = 0, 0 --Fourier coefficients in bin k.
for n = 0, N - 1 do
local cos, sin = Twiddle( n * k )
x = x + seq[n + 1] * cos
y = y + seq[n + 1] * sin
end
dre[k + 1], dim[k + 1] = x / N , y / N
end
end
end
--Minimal-oscillation interpolant of a real function from its discrete Fourier coefficients.
Interpolate = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = re[1]
for k = 1, N / 2 do
local c, s = math.cos( x * k ), math.sin( x * k )
y = y
+ c * re[k + 1]
- s * im[k + 1]
+ c * re[N - k + 1]
+ s * im[N - k + 1]
end
return y
end
--Derivative of the interpolation.
Derivative = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = 0
for k = 1, N / 2 do
local c, s = k * math.cos( x * k ), k * math.sin( x * k )
y = y
- c * im[k + 1]
- s * re[k + 1]
+ c * im[N - k + 1]
- s * re[N - k + 1]
end
return y
end
--Second derivative of the interpolation.
SecondDerivative = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = 0
for k = 1, N / 2 do
local c, s = -k * k * math.cos( x * k ), -k * k * math.sin( x * k )
y = y
+ c * re[k + 1]
- s * im[k + 1]
+ c * re[N - k + 1]
+ s * im[N - k + 1]
end
return y
end
--Bandlimited impulse located at angle theta.
local function AliasedSinc( theta, x )
local n, d = math.sin( N * 0.5 * ( x - theta ) ), N * math.sin( 0.5 * ( x - theta ) )
if d < 0.0001 and d > -0.0001 then return 1.0 end
return n / d
end
--Apply bandlimited impulse to wave, adjust free parameters according to game state.
local function OnImpact( impact, level )
IMPULSESIZE = 10.0 * ( level - 2.0) / 120.0
SOUNDSPEED = 25 - 10 * level / 120
DAMPING = 0.01 * ( 1.0 - 0.4 * level / 120 )
--Apply bandlimited impulse.
local r = cur.radii
local theta = impact.th
local magnitude = IMPULSESIZE * impact.speed
local dt = 1.0 / 120.0
for i = 0, N - 1 do
r[ i + 1 ] = r[ i + 1 ] + dt * magnitude * AliasedSinc( theta, 2.0 * math.pi * i / N )
end
--We've updated the positions, now we need to take a DFT
--in order to get the bandlimited second spatial derivative.
cur:DFT()
return Integrate( dt )
end
local mt = { __index = {
DFT = DFT,
Interpolate = Interpolate,
Derivative = Derivative,
SecondDerivative = SecondDerivative,
Impulse = Impulse } }
local function Wave( )
local t = {
--radii[k] = radius of point on curve at angle (k - 1) / N
radii = {},
--SPACE DFT of radius function (which is periodic)
dftre = {},
dftim = {},
}
for i = 1, N do
t.radii[i] = 1.0
end
DFT( t )
return setmetatable(t, mt)
end
local function Draw( score )
-- Blue circle.
love.graphics.setColor( 91 / 255, 206 / 255, 250 / 255 )
shader:send( "re", unpack( cur.dftre ) )
shader:send( "im", unpack( cur.dftim ) )
shader:send( "score", score )
love.graphics.setShader( shader )
love.graphics.circle("fill", 0, 0, 1.5)
local t = love.timer.getTime()
love.graphics.setShader( )
-- Debug dots.
--[[
for i = 1, N do
local th = ( i - 1 ) * 2.0 * math.pi / N
local cx, cy = cur.radii[i] * math.cos( th ), cur.radii[i] * math.sin( th )
love.graphics.setCanvas()
love.graphics.setColor( 0, 0, 0, 0.5 )
love.graphics.circle( "fill", cx, cy, 0.02 )
for k = 0.1, 1.0, 0.1 do
--Interpolant.
love.graphics.setColor( 1.0, 0, 0, 0.7 )
th = ( i - 1 + k ) * 2.0 * math.pi / N
local r = cur:Interpolate( th )
local x, y = r * math.cos( th ), r * math.sin( th )
love.graphics.circle( "fill", x, y, 0.01 )
--love.graphics.circle( "fill", th / math.pi - 1.0 , r, 0.01)
--First derivative.
--love.graphics.setColor( 0, 1.0, 0, 0.7 )
r = cur:Derivative( th )
x, y = r * math.cos( th ), r * math.sin( th )
--love.graphics.circle( "fill", x, y, 0.01 )
love.graphics.circle( "fill", th / math.pi - 1.0, r, 0.01)
--Second derivative.
love.graphics.setColor( 0, 0, 1.0, 0.7 )
r = cur:SecondDerivative( th )
x, y = r * math.cos( th ), r * math.sin( th )
--love.graphics.circle( "fill", x, y, 0.01 )
love.graphics.circle( "fill", th / math.pi - 1.0, r, 0.01)
love.graphics.setColor( 1.0, 1.0, 1.0, 0.2 )
love.graphics.circle( "fill", 2.0 * ( i + k ) / N - 1.2, 0, 0.02 )
end
end
love.graphics.setColor( 1, 1, 1, 0.5 )
local r = cur:Interpolate( t )
local x, y = r * math.cos( t ), r * math.sin( t )
love.graphics.circle( "fill", x, y, 0.02 )]]
end
local function Update( )
--Deep copy of current state to old state.
for name, t in pairs( cur ) do
for i = 1, N do
old[name][i] = t[i]
end
end
--Deep copy of new state to current state.
for name, t in pairs( new ) do
for i = 1, N do
cur[name][i] = t[i]
end
end
end
Integrate = function( step )
for i = 1, N do
local rxx = cur:SecondDerivative( math.pi * 2.0 * ( i - 1 ) / N )
local r = ( 1.0 - DAMPING ) * ( 2.0 * cur.radii[i] - old.radii[i] + step * step * SOUNDSPEED * rxx ) --Verlet
+ DAMPING --Damping: oscillate toward 1.
--Avoid explosions.
r = math.max( 0.5, math.min( r, 4.0 ) )
new.radii[i] = r
end
new:DFT( )
end
local function DetectCollision( px, py, vpx, vpy )
local xi, xf = px + vpx * step, py + vpy * step
end
local function Reset()
IMPULSESIZE = 1 / 10.0
SOUNDSPEED = 5.0
DAMPING = 0.1 / 1
old = Wave()
cur = Wave()
new = Wave()
end
Reset()
return {
Reset = Reset,
Update = Update,
Integrate = Integrate,
OnImpact = OnImpact,
DetectCollision = DetectCollision,
Draw = Draw,
Current = Current,
Next = Next,
}