# =============================================================================
# phaseshaper_apps.py -- Phaseshaper applications
# (section 3 of the paper)
#
# 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 and Matplotlib v0.99.1.1
# uses pylab for straight-forward matlab porting
# tab settings: 3
# =============================================================================
from scipy.signal import *
import matplotlib
import pylab

from phaseshapers import *		# phaseshaper implementations


# -----------------------------------------------------------------------------
# Synthesis
# phaseshaper implementations in phaseshapers.py
# -----------------------------------------------------------------------------

fs = 44100.							# sample rate
t = arange(0,fs)					# time vector

# -- initial phase signals
# -- f0 = 441,392,1245 Hz (one signal each)
phi1 = zeros(len(t))		# phase
phi2 = zeros(len(t))
phi3 = zeros(len(t))
si1 = 441./fs				# phase increments
si2 = 392./fs
si3 = 1245./fs
P1 = int(1/si1)			# periods
P2 = int(1/si2)
P3 = int(1/si3)
phi1[0] = phi2[0] = phi3[0] = 0.0
for n in range(1,len(t)):
	phi1[n] = (phi1[n-1] + si1) % 1
	phi2[n] = (phi2[n-1] + si2) % 1
	phi3[n] = (phi3[n-1] + si3) % 1

# -- waveslice (sect 3.1 and 3.2)
slice_phase = g_lin(phi3,a1=0.25)
slice_triv = g_b(sin(2*pi*slice_phase))	# trivial
slice_aa = g_mods(slice_triv,si3,-2)		# soft modulo (alias-suppressed)

# -- oscillator sync (sect 3.3)
hardsync = g_b(g_ramp(phi1,a1=2.5))
softphase = g_tri(phi1,a1=1.25)
softsync = g_b(s_tri(softphase))

# -- triangle modulation (sect 3.5)
atm = 0.82
trimod_phase = atm*g_b(g_tri(phi1))
trimod = 2*(trimod_phase-ceil(trimod_phase-0.5))

# -- supersaw (sect 3.6, three snapshots capturing the LFO-modulated evolution)
m2 = 0.88
a1 = 1.5
xs = g_lin(phi1,a1)
supersaw_phase_a = modm(xs,0.75) + modm(xs,m2)		# snapshot 1
supersaw_phase_b = modm(xs,0.50) + modm(xs,m2)		# snapshot 2
supersaw_phase_c = modm(xs,0.25) + modm(xs,m2)		# snapshot 3
supersaw_a = g_b(sin(supersaw_phase_a))				# snapshot 1
supersaw_b = g_b(sin(supersaw_phase_b))				# snapshot 2
supersaw_c = g_b(sin(supersaw_phase_c))				# snapshot 3

# -- half/full-cycle sinusoids (sect 3.7)
# var-slope
w1 = 0.50
w2 = 0.11	# for causal g_pulse
pulse1 = g_pulse(phi2,w1,P2)
vslope1 = 0.5*phi2*(1-pulse1)/w1 + pulse1*(phi2-w1)/(1-w1)
sin1 = sin(2*pi*vslope1)
pulse2 = g_pulse(phi2,w2,P2)
vslope2 = 0.5*phi2*(1-pulse2)/w2 + pulse2*(phi2-w2)/(1-w2)
sin2 = sin(2*pi*vslope2)
# var-triangle
a1 = 1.5
w3 = 0.75
vtri = g_vtri(phi2,w3,a1,0,P2)
sin3 = sin(2*pi*vtri)


# -----------------------------------------------------------------------------
# Spectrum (portions by V.Lazzarini)
# -----------------------------------------------------------------------------

N = 16384
win = chebwin(N, 100)
scal = N*sqrt(mean(win**2))
bins = t[0:N/2]*fs/N

def calcSpectrum(sig):
	spec = fft(win*sig[0:N])
	mags = sqrt(spec[0:N/2].real**2 + spec[0:N/2].imag**2)
	norm = 20*log10(mags/scal)
	norm = norm - max(norm)
	return norm


# -----------------------------------------------------------------------------
# Plotting
# -----------------------------------------------------------------------------

def preplot(id,signal=True):
	p = pylab.subplot(id)
	p.grid(True)
	if signal:
		ticks = matplotlib.ticker.FixedLocator([-1,-0.5,0,0.5,1])
		p.yaxis.set_major_locator(ticks)
	return p

def postplot(p,signal=True):
	if signal:
		pylab.ylim(-1.1, 1.1)
		pylab.xlim(0,150)
		pylab.ylabel('Level')
	else:
		pylab.ylim(-80,0)
		pylab.xlim(0,fs/2)
		ticks = matplotlib.ticker.FixedLocator([-80,-60,-40,-20,0])
		p.yaxis.set_major_locator(ticks)
		p.xaxis.set_ticklabels(['','5','10','15','20'])
		pylab.ylabel('Magnitude (dB)')
		pylab.xlabel('Frequency (kHz)')


# -- waveslices (sect 3.1: trivial)
fig = pylab.figure(figsize=(7,9))
fig.canvas.set_window_title('Waveslices')
p11 = preplot(411)
p11.plot(slice_phase, "r", lw=1)
p11.plot(slice_triv, "k", lw=2)
postplot(p11)
pylab.xlabel('Time (samples, f0=1245 Hz, fs=44.1 kHz)')
p12 = preplot(412,False)
p12.plot(bins, calcSpectrum(slice_triv), "k", lw=1)
postplot(p12,False)

# -- waveslices (sect. 3.2: aliasing-suppressed)
p13 = preplot(413)
p13.plot(slice_phase, "r", lw=1)
p13.plot(slice_aa, "k", lw=2)
postplot(p13)
pylab.xlabel('Time (samples, f0=1245 Hz, fs=44.1 kHz)')
p14 = preplot(414,False)
p14.plot(bins, calcSpectrum(slice_aa), "k")
postplot(p14,False)
pylab.subplots_adjust(hspace=0.37)

# -- oscillator sync (sect 3.3)
# hardsync
fig = pylab.figure(figsize=(6,4.5))
fig.canvas.set_window_title('Oscillator Sync + Triangle Modulation')
p21 = preplot(311)
p21.plot(hardsync, "k", lw=2)
postplot(p21)
pylab.xlim(0,301)
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')
# softsync
p22 = preplot(312)
p22.plot(softphase[50:], "r", lw=1)
p22.plot(softsync[50:], "k", lw=2)
postplot(p22)
pylab.xlim(0,301)
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')

# -- triangle modulation (sect 3.5)
p41 = preplot(313)
p41.plot(trimod_phase[74:], "r", lw=1)
p41.plot(trimod[74:], "k", lw=2)
postplot(p41)
pylab.xlim(0,301)
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')
pylab.subplots_adjust(bottom=0.13,hspace=0.5)

# -- supersaw (sect 3.6)
fig = pylab.figure(figsize=(6,6.5))
fig.canvas.set_window_title('Supersaw')
# snapshot 1
p51 = preplot(311)
p51.plot(supersaw_phase_a, "r", lw=1)
p51.plot(supersaw_a, "k", lw=2)
postplot(p51)
pylab.ylim(-1.6,1.6)
pylab.xlim(0,301)
p51.yaxis.set_major_locator(matplotlib.ticker.FixedLocator([-1.5,-1,-0.5,0,0.5,1,1.5]))
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')
# snapshot 2
p52 = preplot(312)
p52.plot(supersaw_phase_b, "r", lw=1)
p52.plot(supersaw_b, "k", lw=2)
postplot(p52)
pylab.ylim(-1.6,1.6)
pylab.xlim(0,301)
p52.yaxis.set_major_locator(matplotlib.ticker.FixedLocator([-1.5,-1,-0.5,0,0.5,1,1.5]))
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')
# snapshot 3
p53 = preplot(313)
p53.plot(supersaw_phase_c, "r", lw=1)
p53.plot(supersaw_c, "k", lw=2)
postplot(p53)
pylab.ylim(-1.6,1.6)
pylab.xlim(0,301)
p53.yaxis.set_major_locator(matplotlib.ticker.FixedLocator([-1.5,-1,-0.5,0,0.5,1,1.5]))
pylab.xlabel('Time (samples, f0=441 Hz, fs=44.1 kHz)')
pylab.subplots_adjust(hspace=0.4)

# -- full and half-cycle sinusoids (sect 3.7)
fig = pylab.figure(figsize=(13,6))
fig.canvas.set_window_title('Half/Full-cycle Sinusoids')
# fig. 14a
p61 = preplot(321)
p61.plot(vslope1, "r", lw=1)
p61.plot(sin1, "k", lw=2)
postplot(p61)
pylab.xlim(0,3*P2)
pylab.xlabel('Time (samples, f0=392 Hz, fs=44.1 kHz)')
p62 = preplot(322,False)
p62.plot(bins, calcSpectrum(sin1), "k")
postplot(p62,False)
# fig. 14b
p63 = preplot(323)
p63.plot(vslope2, "r", lw=1)
p63.plot(sin2, "k", lw=2)
postplot(p63)
pylab.xlim(0,3*P2)
pylab.xlabel('Time (samples, f0=392 Hz, fs=44.1 kHz)')
p64 = preplot(324,False)
p64.plot(bins, calcSpectrum(sin2), "k")
postplot(p64,False)
# fig. 15
p65 = preplot(325)
d = (int(w3*P2))-1
p65.plot(vtri[d:], "r", lw=1)
p65.plot(sin3[d:], "k", lw=2)
postplot(p65)
pylab.xlim(0,3*P2)
pylab.xlabel('Time (samples, f0=392 Hz, fs=44.1 kHz)')
p66 = preplot(326,False)
p66.plot(bins, calcSpectrum(sin3), "k")
postplot(p66,False)
pylab.subplots_adjust(hspace=0.4)

pylab.show()