Phaseshaping Osc Algorithms

[jamshark70]

As an exercise, and also to illustrate some points about efficient SynthDef writing, here is a revision to Maths2 that also fixes the bug about scalar linExp values.

Edit: I’m being picky about it for a couple of reasons – mainly to demonstrate some things that we can think about while writing synths. Sam explained later that this was written very quickly on low battery without access to power, so it wasn’t finished. That’s fine! I hope it’s helpful that I went ahead with some updates.

The main points:

• Multiplication is faster than division: so, x * SampleDur.ir is more efficient than x / SampleRate.ir.
• SC does not automatically fold duplicate UGens down into one. If you write SampleDur.ir 10 times, there will be 10 instances of SampleDur.
  • This is even more important for operators. If you write linExp > 0.5 six times, it will actually do the identical comparison operation six times. It’s much more efficient to save the result in a variable and reuse the variable in multiple places.
• The linear interpolation formula (y - x) * interp + x has one fewer multiplication than (y * interp) + (x * (1 - interp)).
  • In general, it’s worth looking for algebraic reductions: e.g., x.neg > 0 is the same as x < 0, except that the latter renders into a single UGen where the former requires two.

Maths2 {
	*ar { |rise = 0.1, fall = 0.1, linExp = 0.5, loop = 1, plugged = 0, trig = 0|

		var freq = (rise.clip(0.001, 10*60) + fall.clip(0.001, 10*60)).reciprocal;
		var width = rise.clip(0.001, 10*60) / (rise.clip(0.001, 10*60) + fall.clip(0.001, 10*60));

		// avoid multiple instances of the same UGen,
		// if it's known that they're all the same.
		// also dividing by SampleRate is slower than multiplying by SampleDur
		var sampleDur = SampleDur.ir;

		var plugTrig = Trig1.ar(1 - plugged, sampleDur);
		var loopTrig = Trig1.ar(loop, sampleDur);

		var eof, eor;
		var maths, maths2, interp, isExp;

		var phasor = Phasor.ar(Silent.ar + loopTrig + plugTrig, 2 * freq * sampleDur, -1, 1, -1);

		var phasorTrig, latchTrig, postEnv;
		var inTrig, phasor2, postEnv2;

		var pluggedIsVariable = plugged.rate != \scalar;
		var needMaths = pluggedIsVariable or: { plugged < 1 };
		var needMaths2 = pluggedIsVariable or: { plugged >= 1 };

		if(needMaths) {
			phasorTrig = Trig1.ar(0.5 - phasor, sampleDur) + EnvGen.ar(Env([0, 0, 1, 0], [sampleDur, sampleDur, sampleDur]));
			latchTrig = (phasorTrig + (DelayN.ar(loopTrig, 0.01, 0.01))).clip(0, 1);
			postEnv = (Latch.ar(K2A.ar(loop), latchTrig) > 0);

			phasor = phasor.linlin(-1, 1, width.neg, 1-width);
			maths = phasor.bilin(0, width.neg, 1-width, 0, -1, 1);
			maths = 1 - maths.abs;
			maths = maths * postEnv;
		};

		if(needMaths2) {
			inTrig = Trig1.ar(trig, sampleDur);
			phasor2 = Phasor.ar(inTrig, 2 * freq * sampleDur, -1, 1, -1);
			postEnv2 = SetResetFF.ar(Delay1.ar(inTrig), Trig1.ar(0.5 - phasor2, sampleDur));
			phasor2 = phasor2.linlin(-1, 1, width.neg, 1 - width);
			maths2 = phasor2.bilin(0, width.neg, 1 - width, 0, -1, 1);
			maths2 = 1 - maths2.abs;
			maths2 = maths2 * postEnv2;
		};

		if(pluggedIsVariable) {
			maths = Select.ar(plugged, [maths, maths2]);
		} {
			if(plugged >= 1) { maths = maths2 };
		};

		isExp = linExp > 0.5;  // don't repeat this 3 times

		// could reduce here as well, if linExp is scalar
		// but I'm tired now ;-p and these aren't expensive
		interp = Select.kr(isExp, [linExp.linlin(0, 0.5, 1, 0), linExp.linlin(0.5, 1, 0, 1)]);
		maths = Select.ar(isExp, [maths - 1, maths]);
		// maths = (maths ** 8 * interp) + (maths * (1 - interp));
		// - * + is more efficient than *, - * +
		maths = ((maths ** 8) - maths) * interp + maths;
		maths = Select.ar(isExp, [maths + 1, maths]);

		if(pluggedIsVariable) {
			// again, repeated multiplications...
			postEnv = phasor * postEnv;
			postEnv2 = phasor2 * postEnv2;
			// x.neg > 0, same as x < 0, save one unary op
			eof = Select.ar(plugged, [postEnv < 0, postEnv2 < 0]);
			eor = Select.ar(plugged, [postEnv > 0, postEnv2 > 0]);
		} {
			if(plugged < 1) {
				postEnv = phasor * postEnv;
				eof = postEnv < 0;
				eor = postEnv > 0;
			} {
				postEnv2 = phasor2 * postEnv2;
				eof = postEnv2 < 0;
				eor = postEnv2 > 0;
			};
		};

		^[maths, eof, eor]
	}
}

hjh

[dietcv]

thanks alot for all the input @jamshark70 @Sam_Pluta. i was trying to get the tilted triangle out of maths2. looked at the source code for the pseudo uGen. but couldnt figure it out. any suggestions?

(
{
	var linExp = \linExp.kr(0.5);
	var rise = \rise.kr(0.01);
	var fall = \fall.kr(0.01);

	var freq = \freq.kr(150);
	var width = \width.kr(0.25);
	var trig = Impulse.ar(freq);
	var bufFrames = BufFrames.ir(b);
	var sampleDur = SampleDur.ir;

	var maths = Maths2.ar(
		rise: rise,
		fall: fall,
		linExp: linExp,
		loop: 0,
		plugged: 1,
		trig: trig
	)[0];

	var tilted_tri_maths = maths.wrap(0,1);

	var phasor = Phasor.ar(trig, freq * sampleDur, width.neg, 1-width, width.neg);
	var tilted_tri = (phasor.bilin(0, width.neg, 1-width, 0, -1, 1).abs*1.5).wrap(0,1);

	[maths, tilted_tri_maths, phasor, tilted_tri]
}.plot
)

[Sam_Pluta]

Notice what I was doing with the tilted tri. I was multiplying the signal by 1.5 before wrapping, so you need to do that here as well. It can’t wrap if it doesn’t go beyond the wrap points, and by nature, it only goes between 0 and 1.

You can get some really cool looking waveforms with this. How they sound is another can of worms!

(
{
	var linExp = \linExp.kr(0.9);
	var rise = \rise.kr(0.001);
	var fall = \fall.kr(0.001);

	var freq = \freq.kr(150);
	var width = \width.kr(0.25);
	var trig = Impulse.ar(freq);
	//var bufFrames = BufFrames.ir(b);
	var sampleDur = SampleDur.ir;

	var maths = Maths2.ar(
		rise: rise,
		fall: fall,
		linExp: linExp,
		loop: 1,
		plugged: 0/*,
		trig: trig*/
	)[0];
	
	var tilted_tri_maths = (maths*1.5).wrap(0,1);

	var phasor = Phasor.ar(trig, freq * sampleDur, width.neg, 1-width, width.neg);
	var tilted_tri = (phasor.bilin(0, width.neg, 1-width, 0, -1, 1).abs*1.5).wrap(0,1);

	[maths, tilted_tri_maths, phasor, tilted_tri]
}.plot
)

Actually, this kind of manipulation (multiply by 1.5, subtract 0.5 and do an absolute value) might sound better in both cases:

(
{
	var linExp = \linExp.kr(0.8);
	var rise = \rise.kr(0.001);
	var fall = \fall.kr(0.001);

	var freq = \freq.kr(150);
	var width = \width.kr(0.25);
	var trig = Impulse.ar(freq);
	//var bufFrames = BufFrames.ir(b);
	var sampleDur = SampleDur.ir;

	var maths = Maths2.ar(
		rise: rise,
		fall: fall,
		linExp: linExp,
		loop: 1,
		plugged: 0/*,
		trig: trig*/
	)[0];
	
	var tilted_tri_maths = (maths*1.5-0.5).abs;

	var phasor = Phasor.ar(trig, freq * sampleDur, width.neg, 1-width, width.neg);
	var tilted_tri = (phasor.bilin(0, width.neg, 1-width, 0, -1, 1).abs*1.5-0.5).abs;

	[maths, tilted_tri_maths, phasor, tilted_tri]
}.plot
)

[dietcv]

thanks a lot. im not sure how to modulate the width in the best possible way for the phase shaping because rise and fall are binded together. but i think it sounds quite nice already

EDIT
In this video NI Reaktor - Vector Phaseshaping Synthesis - YouTube following https://www.researchgate.net/publication/236616428_Vector_Phaseshaping_Synthesis he has one controlknob d for the point of inflection with values between 0-1 and he is building two macros for
1.) phase x < d and 2.) phase x > d
and has two functions for the first macro → 0.5 * x / d and for the second → 0.5 * [1+(x - d / 1 - d) ] (phase distortion). in addition he is replacing the fixed value 0.5 with the variable v for vector phase shaping. i would like to apply this procedure to our tilted tri / slope ramp to build the vector phase shaping osc. so i think i have to change (Select.ar(tri > 0, [0.5 + tri, tri]) right? d is represented in our original code by the variable width right?

is this correct for pd? im confused with .bilin

(
{

	var freq = \freq.kr(150);
	var width = \width.kr(0.10);
	var trig = Impulse.ar(freq);
	var sampleDur = SampleDur.ir;

	var phasor = Phasor.ar(trig, freq * sampleDur, width.neg, 1-width, width.neg);
	//var tri = (phasor.bilin(0, width.neg, 1-width, 0, -1, 1).abs*1.5-0.5).abs;
	var saw = phasor.bilin(0, width.neg, 1-width, 0, -0.5, 1);
	
	var pd = (Select.ar(saw<0, [0.5 * (saw/width), 0.5 * [1+(saw-width/1-width)]])*2pi).sin;

	[saw, pd]
}.plot
)

introducing variable vector for vector phase shaping

(
{

	var freq = \freq.kr(50);
	var vector = MouseY.kr(0.01,3);
	var width = MouseX.kr(0.01,0.99);
	var trig = Impulse.ar(freq);
	var sampleDur = SampleDur.ir;

	var phasor = Phasor.ar(trig, freq * sampleDur, width.neg, 1-width, width.neg);
	var tri = (phasor.bilin(0, width.neg, 1-width, 0, -1, 1).abs*1.5-0.5).abs;
	//var saw = phasor.bilin(0, width.neg, 1-width, 0, -0.5, 1);
	
	//var vps = (Select.ar(saw<0, [vector * (saw/width), (1 - vector) * (saw-width/1-width) + vector ]) * 2pi).sin;
	var vps = (Select.ar(tri<0, [vector * (tri/width), (1 - vector) * (tri-width/1-width) + vector]) * 2pi).sin;

	vps = LeakDC.ar(vps);
	vps!2 * 0.1
}.scope
)

[dietcv]

Im wondering if this is right, it looks exactly like in the paper Figure 2 “VPS waveforms and spectra”:

(
{
	var freq = \freq.kr(500);
	var width = 0.5;
	var vector = 0.85;
	var phasor, vps;
	
	phasor = LFSaw.ar(freq,1).linlin(-1, 1, width.neg, 1-width);
	phasor = phasor.bilin(0, width.neg, 1-width, vector, 0, 1);
	
	vps = (phasor * 2pi).cos.neg;
	
	[phasor, vps]
}.plot;
)

[Sam_Pluta]

Looks good to me. Sounds good to me. Isn’t playing with waveforms fun?

Sam

[dietcv]

haha thanks a lot. yeah its a lot of fun and also a bit frustrating because im missing some of the math. but im learning a lot
I would also like to add the alising suppression described in chapter 3.2 and found this pseudo Ugen classes in Glitch free vowel synthesis test - #2 by Geoffroy. already another difficult task. lets see.

Nearest_Even {
	*ar {
		arg val;
		var val_floor, val_ceil, res, distance;
		val_floor = val.floor;
		val_ceil = val.ceil;
		res = Select.ar (val % 2,
			[ val_floor, val_ceil ],
		);
		distance = (val - res).abs;
		^ [ res, distance ];
	}
}

Nearest_Odd {
	*ar {
		arg val;
		var val_floor, val_ceil, res, distance;
		val_floor = val.floor;
		val_ceil = val.ceil;
		res = Select.ar (val + 1 % 2,
			[ val_floor, val_ceil ],
		);
		distance = (val - res).abs;
		^ [ res, distance ];
	}
}

Crossfade_Formant {
	*ar {
		arg phasor = 0, phase_mod = 0, harm = 1, pm_index = 1, amp = 1; // lag = 0;
		var harm_even, harm_odd, sig_even, sig_odd, sig, phase;

		phase = phase_mod * pm_index;

		harm_even = Nearest_Even.ar (harm);
		sig_even = cos (phasor * 2pi * harm_even[0] + phase);

		harm_odd = Nearest_Odd.ar (harm);
		sig_odd = cos (phasor * 2pi * harm_odd[0] + phase);

		sig = XFade2.ar (sig_even, sig_odd, harm_even[1] * 2 - 1) * amp;
		^ sig;
	}
}

[dietcv]

I made a first approach. i think changing the values for the vector is kind of smooth now, but im not sure if its the exact same approach presented in the paper. any further ideas?

(
~nearest_even = {|val|
	var val_floor, val_ceil, res, distance;
	val_floor = val.floor;
	val_ceil = val.ceil;
	res = Select.ar (val % 2,
		[ val_floor, val_ceil ],
	);
	distance = (val - res).abs;
	[res, distance];
};

~nearest_odd = {|val|
	var val_floor, val_ceil, res, distance;
	val_floor = val.floor;
	val_ceil = val.ceil;
	res = Select.ar (val + 1 % 2,
		[ val_floor, val_ceil ],
	);
	distance = (val - res).abs;
	[res, distance];
};
)

(
SynthDef(\vector, {
	arg out=0, amp=0.25, freq=150;
	var horizontal = MouseX.kr(0.01,0.99);
	var vertical = K2A.ar(MouseY.kr(1,10)); 
	var sig, sig_even, sig_odd, phasor, phasor_even, phasor_odd, harm_even, harm_odd;

	harm_even = ~nearest_even.(vertical);
	harm_odd = ~nearest_odd.(vertical);
	
	phasor = LFSaw.ar(freq,1).linlin(-1, 1, horizontal.neg, 1-horizontal);
	phasor_even = phasor.bilin(0, horizontal.neg, 1-horizontal, harm_even[0], 0, 1);
	phasor_odd = phasor.bilin(0, horizontal.neg, 1-horizontal, harm_odd[0], 0, 1);
	
	sig_even = (phasor_even * 2pi).cos.neg;
	sig_odd = (phasor_odd * 2pi).cos.neg;
	
	sig = XFade2.ar (sig_even, sig_odd, harm_even[1] * 2 - 1) * amp;

	sig = Splay.ar(sig);
	Out.ar(out, sig);
}).add;
)

(instrument: \vector, midinote: 31).play;

[dietcv]

I dont know how many hours i already spent on this one, and still not being able to figure it out.
1.)
I think my first approach has somehow the right ingredients but is not yielding the desired result. The Crossfade between even/odd is following somehow the idea presented at the end of chapter 3.1 Synthesising Formants. but the plots dont look right and the sound has a totally different quality to it. Any Ideas?
2.)
Im was also trying to come up with a solution for 3.2 Aliasing Suppression, which should modify the phaseshaper when the phase is inside an incomplete period and should render it as a smooth full-cycle sinusoid like in Fig. 5b.

I think the modified phasor looks like in the paper but the signal is not rendered as a smooth full-cycle sinusoid at the inflection point. Any ideas?

(
{
	arg freq = 392;
	var horizontal = \horizontal.kr(0.80);
	var vertical = \vertical.ar(2.2);
	var cos, phasor;

	phasor = LFSaw.ar(freq/2,1).range(horizontal.neg, 1-horizontal);
	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);
	phasor = Select.ar(phasor > vertical.floor, [phasor, phasor.abs.wrap(0,1)]);

	cos = (phasor * 2pi).cos.neg;

	[phasor, cos]
}.plot
)

[Sam_Pluta]

Not sure if this is “right” but just adding an LPF fixes the look of it:

(
{
	arg freq = 392;
	var horizontal = \horizontal.kr(0.80);
	var vertical = \vertical.ar(2.2);
	var cos, phasor;

	phasor = LFSaw.ar(freq/2,1).range(horizontal.neg, 1-horizontal);
	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);
	phasor = Select.ar(phasor > vertical.floor, [phasor, phasor.abs.wrap(0,1)]);

	cos = LPF.ar((phasor * 2pi).cos.neg, (freq*16).clip(0,SampleRate.ir/2));
	//cos = (phasor * 2pi).cos.neg;
	
	[phasor, cos]
}.plot
)

[dietcv]

hey thanks for your help again.

running your code i get:

ERROR: Primitive ‘_ClipInt’ failed.

[Sam_Pluta]

Huh. I don’t get that. Nor have I ever seen that. Try this:

cos = LPF.ar((phasor * 2pi).cos.neg, (freq*16).clip(20.0,SampleRate.ir/2.0));

[dietcv]

waoh thats interesting i get still the same error. ERROR: Primitive ‘_ClipInt’ failed.

[Sam_Pluta]

cos = LPF.ar((phasor * 2pi).cos.neg, min((freq*16.0), SampleRate.ir/2));

??

[dietcv]

this is not giving an error. thanks for the hack but I really dont understand why phasor looks the same but does not have any effect on cos. just one missing piece i guess… these two equations (11) and (12) are giving me sleepless nights haha

[Sam_Pluta]

This is because it is basically just using the first part of the cos curve and the last part of the cos curve - the beginning of the upward slope and the end of the downward slope. It isn’t plotting the full curve. You are just getting the circled part of the drawing:

[Sam_Pluta]

OK. It gets weird if the vertical isn’t between 2 and 3, but this is the jam. The thing that was probably messing you up is the smaller cosine. That part of the signal was basically (after the math) going from 0 to 0.2. But you need it to go from (0 to 1)*pi in order for it to give you a full cosine, thus the division by vertical-2. Then, after getting the cosine, you need to scale the signal back down by multiplying by vertical-2.

(
{
	arg freq = 300;
	var horizontal = \horizontal.kr(0.80);
	var vertical = \vertical.ar(2.2);
	var cos, phasor;

	phasor = LFSaw.ar(freq/2,1).range(horizontal.neg, 1-horizontal);
	
	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);
	cos = Select.ar(phasor > vertical.floor, 
		[
			(phasor * 2pi).cos.neg, 
			((phasor.abs.wrap(0,1)/(vertical-2)*pi).cos*(vertical-2)).neg-(3-vertical)
		]
	);
	
	[phasor, cos+1]
}.plot
)

[dietcv]

thats really great i cant thank you enough its sounding really nice.

some additonal things:
1.) the .abs before .wrap is not necessary right?
2.) should I always add 1 to cos for the final output?
3.) should freq always be devided by 2? has this something to do with the nyquist frequency? in other examples where a Phasor was used for ringsync for example ive seen that freq was multiplied by 2 and the fundamental was just freq. this is confusing me a bit.
3.) considering the fact that vertical should be always between 2-3 it will be kind of tricky to implement multiple inflection to get square like shapes without aliasing described in chapter 3.4 right? Im also not sure how this can be done in general.

thanks so much

EDIT: I also found this thread where others tried to rebuild the alias suppression in reaktor Vector phase shaping alias suppression? | NI Community Forum and the original sound examples + python code (unfortunately i dont know anything about phython) from the authors VPS - Vector Phaseshaping Synthesis

[dietcv]

It gets weird if the vertical isn’t between 2 and 3, but this is the jam

waoh, i found out when you change the values for the scaling vertical-2 to vertical-3 and 3-vertical to 4-vertical, you can go to values between 3 and 4 instead of 2 and 3 for the vertical and so on indefinitely. is there a way to have the scaling behave dynamically related to the values for vertical? Then it would be possible to do a smooth sweep from vertical 0.5 - n i guess or? and introduce some more formants.

(
{
	arg freq = 300;
	var horizontal = \horizontal.kr(0.50);
	var vertical = \vertical.ar(3.5);
	var cos, phasor;

	phasor = LFSaw.ar(freq/2,1).range(horizontal.neg, 1-horizontal);
	
	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);
	cos = Select.ar(phasor > vertical.floor, 
		[
			(phasor * 2pi).cos.neg, 
			((phasor.abs.wrap(0,1)/(vertical-3)*pi).cos*(vertical-3)).neg-(4-vertical)
		]
	);
	
	[phasor, cos+1]
}.plot
)

[dietcv]

EDIT: The code now looks perfect on the plot for every vertical value.
But the transitions are really rough, there is a discontinuity for every integer value because of the case vertical == vertical.floor and its division trough zero. so ive implemented another Select.ar with BinaryOpUGen('==', vertical, vertical.floor. i think when this condition is met the signal should be rendered as a normal cosine, because no anti-aliasing is needed.
but this is probably wrong, because now the signal is not rendered as a smooth full-cycle sinusoid when the phase is inside an incomplete period any more. any ideas?

(
{
	arg freq = 392;
	var horizontal = \horizontal.kr(0.80);
	var vertical = \vertical.kr(2.2);
	var cos, phasor;

	phasor = LFSaw.ar(freq/2,1).range(horizontal.neg, 1-horizontal);

	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);
	cos = Select.ar(phasor > vertical.floor,
		[
			(phasor * 2pi).cos.neg,
			Select.ar(BinaryOpUGen('==', vertical, vertical.floor),
				[
					(phasor * 2pi).cos.neg,
					((phasor.wrap(0,1)/(vertical-vertical.floor)*pi).cos*(vertical-vertical.floor)).neg-(vertical.ceil-vertical)
				]
			);
		]
	);
}.plot
)

EDIT: i tried the crossfade attempt once more because i couldnt figure out how to solve the case vertical = vertical.floor, the result is smooth for sweeping the vertical value but i dont know if its correct and it probably just fixes the vertical value, the aliasing is still prominent for extreme horizontal values.

(
~nearest_even = {|val|
	var val_floor, val_ceil, res, distance;
	val_floor = val.floor;
	val_ceil = val.ceil;
	res = Select.ar (val % 2,
		[ val_floor, val_ceil ],
	);
	distance = (val - res).abs;
	[res, distance];
};

~nearest_odd = {|val|
	var val_floor, val_ceil, res, distance;
	val_floor = val.floor;
	val_ceil = val.ceil;
	res = Select.ar (val + 1 % 2,
		[ val_floor, val_ceil ],
	);
	distance = (val - res).abs;
	[res, distance];
};
)

(
{
	arg out=0, pan=0, amp=0.25, freq=110;

	var horizontal = MouseX.kr(0.01,0.99);
	var vertical = K2A.ar(MouseY.kr(1.0,10.0));

	var vertical_even = ~nearest_even.(vertical);
	var vertical_odd = ~nearest_odd.(vertical);
	var cos, phasor, sig;

	vertical = [vertical_even[0], vertical_odd[0]];

	phasor = Phasor.ar(0, freq/2 * SampleDur.ir, horizontal.neg, 1-horizontal, horizontal.neg);
	phasor = phasor.bilin(0, horizontal.neg, 1-horizontal, vertical, 0, 1);

	cos = (phasor * 2pi).cos.neg;

	sig = XFade2.ar(cos[0], cos[1], vertical_even[1] * 2 - 1);

	sig = LeakDC.ar(sig);
	sig = Pan2.ar(sig, pan, amp);
	Out.ar(out, sig);
}.scope;
)

s.freqscope;

EDIT: the Phython Code from the Authors says something like this:

f0 = 500.
fs = 44100.
w0 = 2*pi*f0/fs
N = int(fs/f0)
t = arange(0,44100)
T = f0/fs
P = 1/T

d = 0.80
v = 2.20

# -- phase distortion function
def vps(x,d,v, suppress):
if x < d: p = (v*x)/d
else: p = (1-v)*(x-d)/(1-d) + v
# aliasing reduction
c1 = c2 = 0
if suppress:
if p > int(v):
c = v % 1
if 0 < c and c <= 0.5:
p = 0.5 * (p % 1) / c
c1 = 0.5 * (1 - cos(2*pi*c))
c2 = c1-1
else:
p = (p % 1) / c
if p > 0.5:
c1 = 0.5 * (1 + cos(2*pi*c))
c2 = 1-c1
return p,c1,c2

# -- signal generation
def process(d,v, suppress):
y = zeros(len(t))
x = zeros(len(t))
phi = 0
for n in range(0,len(t)):
p,c1,c2 = vps(phi,d,v, suppress)
y[n] = -cos(2*pi*p) # waveshaping
x[n] = p
if c1 != 0:
y[n] = y[n]*c1 + c2
phi += T
if phi > 1: phi -= 1
return y,x

# -- spectrum
def spectrum(s):
NFFT = 16384*2
NFFT2 = 16384
NWIN = 16384
win = chebwin(NWIN, 120)
win = append(win, zeros(NFFT-NWIN))
scal = NFFT*sqrt(mean(win**2))
spec = fft(win*s[0:NFFT])
mags = sqrt(spec[0:NFFT2].real**2 + spec[0:NFFT2].imag**2)
norm = 20*log10(mags/scal)
spec = norm - max(norm)
return spec


# -- main
y1,x1 = process(d,v, False)
y2,x2 = process(d,v, True)