# =============================================================================
# phaseshapers.py -- Phaseshaper implementations
#
# J.Kleimola, V.Lazzarini, J.Timoney, V.Valimaki,
# "PHASESHAPING OSCILLATOR ALGORITHMS FOR MUSICAL SOUND SYNTHESIS",
# SMC 2010, Barcelona, Spain, July 21-24, 2010.
#
# written using SciPy v0.7.1
# tab settings: 3
# =============================================================================
from scipy.signal import *


# -----------------------------------------------------------------------------
# Entities
# -----------------------------------------------------------------------------
def mod1(x):	return (x % 1)
def modm(x,m):	return (x % m)
def g_b(x):		return (2*x - 1)
def g_u(x):		return (0.5*x + 0.5)


# -----------------------------------------------------------------------------
# Elementary phaseshapers
# -----------------------------------------------------------------------------

# -- linear transformations
def g_lin(x, a1=1, a0=0):		return (a1*x + a0)
def g_ramp(x, a1=1, a0=0):		return mod1( g_lin(x,a1,a0) )
def g_tri(x, a1=1, a0=0):		return mod1( g_lin(abs(g_b(x)),a1,a0) )

# -- phase-shifted differences
def g_pulse(x, w=0.5, P=100):
	d = int(w*P)
	xdiff = zeros(len(x))
	for n in range(0,len(x)):
		xdiff[n] = x[n] - x[n-d] + (1-w)		# causal
	#  xdiff[n] = x[n] - x[n+d] + w			# non-causal
	return xdiff

def s_vtri(x, w=0.5, P=100):
	d = ceil(w*P+0.5)
	xdiff = zeros(len(x))
	for n in range(d,len(x)):
		x1 = 2*x[n]-1
		xd = 2*x[n-d]-1
		xdiff[n] = x1*x1-xd*xd
	at = 1/(8.*(w-w*w))
	bt = 0.5
	s_vtri = g_lin(xdiff,at,bt)
	s_vtri[0:d] = 0
	return s_vtri

def g_vtri(x, w=0.5, a1=1, a0=0, P=100):
	vtri = s_vtri(x,w,P)
	return mod1(g_lin(vtri,a1,a0))

# -- ripples
def g_ripple(x, m=0.0):
	return x + modm(x,m)


# -----------------------------------------------------------------------------
# Piecewise linear definitions
# -----------------------------------------------------------------------------

# -- triangle
def s_tri(x):
	s = zeros(len(x))
	for n in range(0,len(x)):
		if x[n] < 0.5:	s[n] = 2*x[n]
		else:				s[n] = 2-2*x[n]
	return s


# -----------------------------------------------------------------------------
# Aliasing suppression
# -----------------------------------------------------------------------------

# x = signal
# T = phase increment (f0/fs)
# h = transition height (negative for falling transitions)
# (after Valimaki + Huovilainen, IEEE SPM, Vol.24, No.2, p.121)
def polyBLEP(x, T, h=-1):
	s = zeros(len(x))					# operate on a copy
	for n in range(0,len(x)):
		p = n*T
		p = p - floor(p)				# phase [0,1)
		if p > (1-T):					# -- before transition -----------------------
			t = p - 1					# fractional phase (negative, in samples)
			t /= T						# fractional phase (-1,0)
			c = 0.5*t*t + t + 0.5	# correction polynomial
			c = c * h					# scale with transition height
			s[n] = x[n] + c			# update sample value
		elif p < T:						# -- after transition ------------------------
			t = p							# fractional phase (positive, in samples)
			t /= T						# fractional phase (0,1)
			c = -0.5*t*t + t - 0.5	# correction polynomial
			c = c * h					# scale with transition height
			s[n] = x[n] + c			# update sample value
		else: s[n] = x[n]				# -- not inside transition area: nop ---------
	return s

# -- soft-modulo using polyBLEP
# -- see above for parameters
def g_mods(x, T, h=-1):
	return polyBLEP(x,T,h)