Barberpole effect with time-varying notch-filters

hey, i was trying to follow this paper https://www.dafx.de/paper-archive/2015/DAFx-15_submission_67.pdf together with the matlab code https://napulen.github.io/media/mumt501/barberpole.zip from this work MUMT 501: Barberpole Phasing and Flanging Illusions | Néstor Nápoles López
to create an barberpole effect for an arbitrary input signal using a cascade of time-varying notch-filters and got stuck on the way. any ideas on this? is this even possible in SC in relation to the blocksize? thanks alot :slight_smile:

(
var getCenterFreqs, getCenterAmps;

getCenterFreqs = { |notch, freqIndex, cycleFreq, f0|
	var exp, centerFreq;
	exp = ( cycleFreq * ( notch - 1) + freqIndex - 1) / cycleFreq;
	centerFreq = f0 * (2 ** exp);
	centerFreq;
};

getCenterAmps = { |notch, freqIndex, notchNums, cycleFreq, mindB, maxdB|
	var numerator, theta, loudness, loudnessCurve, centerGain;
	numerator = ( notch - 1) * cycleFreq + freqIndex - 1;
	theta = 2 * pi * numerator / ( notchNums * cycleFreq );
	loudness = ( maxdB - mindB ) * ( 1 - cos(theta) ) / 2;
	loudnessCurve = mindB + loudness;
	centerGain = 10 ** (loudnessCurve / 20);
	centerGain;
};

SynthDef(\barberpole, {

	var cycleRate = 0.1;
	var mindB = -3.dbamp;
	var maxdB = -20.dbamp;

	var f0 = 20;
	var notchNums = 10;
	var notches = (1..notchNums);
	var sig;

	//sig = Saw.ar(100);
	sig = PinkNoise.ar(1);

	notches.do { |notch, notchNum|
		var freqs, amps;
		var cycleFreq, freqIndex, lfo;

		cycleFreq = floor(SampleRate.ir, cycleRate);
		freqIndex = mod(notchNum, cycleFreq);

		//lfo = LFSaw.ar(cycleFreq, 2 * notch / notchNum).linexp(0, 1, 20, 100);

		freqs = getCenterFreqs.(notch, freqIndex, cycleFreq, f0);
		amps = getCenterAmps.(notch, freqIndex, notchNums, cycleFreq, mindB, maxdB);

		sig = BBandStop.ar(sig, freqs /** lfo*/, \bw.kr(1), amps);
	};

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.25));
	Out.ar(\out.kr(0), sig);
}).play;
)

s.freqscope;

I’d suspect one reason why there hasn’t been much activity on this thread is because the problem isn’t well defined.

Typically troubleshooting goes better with an “expected vs actual” or “wanted - tried - got” model of question.

  • “Expected - actual”: I ran this code. I thought it would do like this (describe what you expected), and this is what I got instead.

  • “Wanted - tried - got” is basically the same, in a different order: I wanted a result like x. I tried to do it with code y, but the result was z.

This question is like “read these other papers to understand what I want” (hm, not really helpful for volunteers who might get involved) and “it’s not doing what I want but I haven’t said what’s bothering me about it.” So it’s not clear what direction to go with it.

Block size matters for LocalIn/LocalOut pairs, or InFeedback. You have a feed-forward chain of filters so block size shouldn’t be an issue.

hjh

Just from skimming the paper, I think the dual flanger and SSB methods would be easier to implement in sc… less filter math to deal with :smiley: Especially the SSB method, sc already has a UGen for SSB frequency shifting (FreqShift) and then you’d just have to convert the -0.5 coefficient for the allpass stages into decay times for AllpassC (the equation is in the help file), I think? Also seems pretty tweakable, you could play around with the coefficient to change the spacing of the notches and the frequency shift amount to change the rate and direction of the effect.

1 Like

thank you for poiting out the obvious.
Im sorry for that. it was early in the morning after sitting there the whole night of trying to convert the matlab code into SC (havent dealt with matlab before).

I was quite happy with the result because it was a difficult task for me trying to understand the filter theory and the matlab code (and still is). The code i have written was at least not producing an error and looked already like im on the right way on the freqscope. But i was not able to formulate a clear question at that moment i think. i will try that now after summarizing the effect im trying to achieve.

The barberpole phasing effect is derived from the classic auditory illusion of the shepard-risset glissando with a continuously rising and falling sweep where the sinusoidal components are equally spaced in the logarithmic frequency scale and each component is one octave higher than the previous one. A bell shaped spectral envelope function (raised-cosine envelope) takes care of the fade-in and fade-outs of the harmonic components.

The barberpole phasing effect as it is described in the paper is inverting the idea of the shepard-risset glissando by replacing the spectral peaks by notches so it can be used on any arbitrary input signal.
For that to work you can use a cascade of notch filters. I have been using the BBandStop Ugen for this and used a cascade of M = 10 Notch-Filters with a center frequency for the first filter at the beginning of each cycle of f0 = 20 Hz as it has been suggested in the paper to ensure the last notch-filter in the chain ends as close to the Nyquist limit as possible. The center frequencies of the notches (like in the shepard-risset glissando) should also be placed in octaves along the frequency spectrum. The amount of attenuation caused by each filter is determined by using an inverted version of the raised-cosine envelope.
the magnitude response of the filter cascade should look like this:

To implement the illusion the paper says: one have to create a rising notch sweep where the center frequencies of the filters move up the spectrum, notches approaching the Nyquist limit will gradually disappear. The center frequencies of the different notches and their respective amplitudes have to be calculated for every single cycle and a low repetition rate for a cycle should be chosen up to about 0.3 Hz at maximum (similiar for the shepard-risset glissando) for the auditory illusion to work. Since the one-octave interval between the notches is preserved at every time step, the center frequencies of the filter wil eventually reach twice their initial value. At this point, we say the system has completed one full cycle. If the filter parameters are continuously reset one time step before the system completes a cycle, the illusion of an infinite filter sweep is generated (mainly referencing the paper here).

In the matlab code the difference equation of a notch-filter: y(n) = b0 * xn - b1 * xnm1 + b2 * xnm2 + a1 * ynm1 - a2 * ynm2 is used to calculate the coefficients of the individual notch-filters for every single cycle.

1.) Can i use a cascade of BBandStop for this to work, or do i have to work with the difference equation itself and calculate the filter coefficients for each cycle? probably possible with Fb1?

2.) Are the formulas ive used for calculating each center frequency and each center amplitude correct (according to the paper) to have a magnitude response of the filter cascade like in figure 4? This question goes towards debugging the code for getCenterFreqs and getCenterAmps. How can i plot the magnitude response of my filter cascade, to make sure its correct?

3.) At the moment nothing is sweeping. Its a static system with just 10 notches. How do i implement the frequency sweep with cycleRate; Do i have to implement some kind of LFO or Impulse with the cycleRate?

4.) Ive converted this from the matlab code: cycleFreq = floor(SampleRate.ir, cycleRate); Does this makes sense in SC or can i just use cycleRate?

thanks a lot :slight_smile:

hey thanks for taking some time :slight_smile: ive also looked at the Dual-Flanger approach and the Freqshift. Trying to use the cascade of Notch filter had also something to do with learning more about filter cascades. I think when the system ive started to write can not be working in the way ive been planning to. then i would go over to using one of the other two approaches.

Apologies if that was an irritating thing to say… I had just noticed that the question sat around for a few days without any response and I was thinking about why.

Long post ahead…

I think so. … except: BBandStop will always fully suppress the center frequency, but you want a gain factor. MidEQ may be what you want instead.

Frequencies will definitely be exponentially distributed in that case (good).

I like to test components in isolation. If you’re not certain that the amplitudes are correct, you could take the function out by itself, look at the numbers, try some plots.

Now, when I look more closely at your inputs to the functions, I see a couple of odd things.

cycleFreq = floor(SampleRate.ir, cycleRate);
freqIndex = mod(notchNum, cycleFreq);

floor is a unary op – there is no second argument. So cycleFreq will only be the sample rate. (Maybe you meant trunc but the purpose of this is still not clear to me.) Then notchNum will always be far, far below even the trunc result. Unless you have tens of thousands of frequencies, this mod it isn’t going to do anything.

So drop those. Just pass integer indices into the functions.

var notches = (1..notchNums); and notches.do { |notch, notchNum| { ... } – in the loop, you will get 1 and 0, then 2 and 1, then 3 and 2, etc… so you don’t need the array at all. You could just do notchNums.do and add 1 where needed… but then, it’s often better to do calculations based on 0 rather than 1. So notch-1 and freqIndex-1 are just working around the “start with 1” anyway… simplify further, get rid of those. Then “notch” and “freqIndex” are the same, so, collapse to one.

Then, after eliminating the spurious distinction between notch and freqIndex, the exp in the frequency function becomes…

exp = (cycleFreq * notch + notch) / cycleFreq;

= ((cycleFreq + 1) * notch) / cycleFreq
= notch * ((cycleFreq + 1) / cycleFreq)
= notch * (1 + (1 / cycleFreq))

cycleFreq is very high, so the fraction will be very small. This means, practically, the exponent is almost exactly equal just to notch – so, there’s really not much point. So I would get rid of that too.

Having simplified, then:

~getCenterFreqs = { |notch, f0|
	f0 * (2 ** notch);
};

But now you don’t have any control over the width of the band. That is, the function parameters don’t satisfy the need (and it’s harder to recognize this fact, because of extra parameters being passed in, and a lot of noise in the formulas, which reduces out once you look at it closely).

So, OK, trying a new approach:

(
~getCenterFreqs = { |notch, notchNums, f0, fn|
	notch.linexp(0, notchNums-1, f0, fn)
};
)

n = (0..9);
o = ~getCenterFreqs.(n, 10, 20, 12000);
o.doAdjacentPairs { |a, b| (b/a).postln }; ""  // all the same, good

The amps can be simplified similarly:

(
~getCenterAmps = { |notch, notchNums, mindB, maxdB|
	var numerator, theta, loudness, loudnessCurve, centerGain;
	theta = 2 * pi * notch / notchNums;
	loudness = ( maxdB - mindB ) * ( 1 - cos(theta) ) / 2;
	loudnessCurve = mindB + loudness;
	centerGain = 10 ** (loudnessCurve / 20);
	centerGain;
};
)

p = ~getCenterAmps.(n, 10, -3, -20);
p.plot;

This looks a bit asymmetrical, but, the barber pole will involve adding a fractional offset, so I think it will be OK.

Then…

(
// different version: MidEQ takes db, not amp! So return dB
~getCenterAmps = { |notch, notchNums, mindB, maxdB|
	var numerator, theta, loudness;
	theta = 2 * pi * notch / notchNums;
	loudness = ( maxdB - mindB ) * ( 1 - cos(theta) ) / 2;
	mindB + loudness;
};

a = SynthDef(\barberpole, {
	var cycleRate = 0.1;
	var mindB = -3;
	var maxdB = -20;

	var f0 = 20, fn = 12000;
	var notchNums = 10;
	var sig;
	var sr = SampleRate.ir;
	
	var rq = MouseX.kr(1, 80, 0).reciprocal;
	
	var freqs = (0..notchNums - 1).linexp(0, notchNums-1, f0, fn);

	//sig = Saw.ar(100);
	sig = PinkNoise.ar(1);

	freqs.do { |freq, i|
		var amp = ~getCenterAmps.(i, notchNums, mindB, maxdB);
		sig = MidEQ.ar(sig, freq, rq, amp);
	};

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.25));
	Out.ar(\out.kr(0), sig);
}).play;
)

The barberpole should be an LFO. I’d test with sines first (because it’s hard to hear a frequency that has been notched out). Actually there’s a mistake in the amp handling here, so this probably doesn’t sound like a real Shepard tone.

(
a = SynthDef(\barberpole, {
	var cycleRate = 0.1;
	var mindB = -3;
	var maxdB = -20;

	var f0 = 20, fn = 12000;
	var notchNums = 10;
	var sig;
	var sr = SampleRate.ir;
	
	var rq = MouseX.kr(1, 80, 0).reciprocal;
	
	var cycle = LFSaw.ar(cycleRate, mul: 0.5, add: 0.5);
	var freqs = (cycle + (0..notchNums - 1)).linexp(0, notchNums, f0, fn);

	//sig = Saw.ar(100);
	// sig = PinkNoise.ar(1);
	sig = 0;

	freqs.do { |freq, i|
		var amp = ~getCenterAmps.(i + cycle, notchNums, mindB, maxdB);
		sig = sig + SinOsc.ar(freq, 0, amp.ampdb)
	};

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.25));
	Out.ar(\out.kr(0), sig);
}).play;
)

Then put the notches back in:

(
a = SynthDef(\barberpole, {
	var cycleRate = 0.1;
	var mindB = -3;
	var maxdB = -20;

	var f0 = 20, fn = 12000;
	var notchNums = 10;
	var sig;
	var sr = SampleRate.ir;
	
	var rq = MouseX.kr(1, 80, 0).reciprocal;
	
	var cycle = LFSaw.ar(cycleRate, mul: 0.5, add: 0.5);
	var freqs = (cycle + (0..notchNums - 1)).linexp(0, notchNums, f0, fn);

	sig = PinkNoise.ar(1);

	freqs.do { |freq, i|
		var amp = ~getCenterAmps.(i + cycle, notchNums, mindB, maxdB);
		sig = MidEQ.ar(sig, freq, rq, amp);
	};

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.25));
	Out.ar(\out.kr(0), sig);
}).play;
)

hjh

thank you very much for collapsing all the formulas to their essentials and showing some methods for debugging. this was really instructive. I havent used the doAdjacentPairs before. thats really handy :slight_smile:

it works beautifully on Pinknoise, as the authors stated this implementation works best on input signals with a dense nearly flat magnitude spectrum.

I just took over the floor(SampleRate.ir, cycleRate) expression from the matlab code, probably with a wrong syntax. Actually exactly this one lead to the confusion about Single Samples, Blocksize, etc.

my initial idea was to use a BPF Saw. this isnt really working. probably better to use the Dual Flanger approach for this.

ive tried my best on the second approach: the synchronized dual flanger.
Im using two CombCs. The LFSaws should have a 90° phase difference and I crossfade between them with XFade2 to avoid the hard reset of the LFSaws as it has been suggested in the paper. it still sounds kind of rough.

1.) are there any better flanger designs?
2.) the paper suggests to use a third-order Lagrangian franctional delay filter to accomodate fractional delay lenghts (to make the notches travel smoothly in frequency). i have no idea about that. somebody else?

(
SynthDef(\flangerSync, {

	var rate = \rate.kr(0.1);
	var minDelay = \minDel.kr(0.001);
	var maxDelay = \maxDel.kr(0.01);
	var decayTime = \decayTime.kr(1);
	var sig, sigA, sigB, lfoA, lfoB, combA, combB, delayTime;

	lfoA = LFSaw.ar(rate, 0, (maxDelay * 0.5) - minDelay, (maxDelay * 0.5) + minDelay);
	lfoB = LFSaw.ar(rate, pi, (maxDelay * 0.5) - minDelay, (maxDelay * 0.5) + minDelay);

	sig = Saw.ar(100);

	combA = CombC.ar(sig, maxDelay, lfoA, decayTime);
	combB = CombC.ar(sig, maxDelay, lfoB, decayTime);

	sig = XFade2.ar(combA, combB, LFTri.kr(rate, pi));

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.1));
	//sig = SafetyLimiter.ar(sig);
	Out.ar(\out.kr(0), sig);
}).play;
)

From LFSaw help: “For efficiency reasons [iphase] is a value ranging from 0 to 2.”

So 90° would be 0.5, not pi.

Actually even if iphase were in radians, then pi would be 180° – so 0.5pi would have been a reasonable guess.

FWIW even after almost 20 years, I still routinely open help files for info on input and output ranges when I’m not sure.

But based on the crossfade, should it perhaps be 180°?

For the crossfade, I’d check to be very sure it’s truly starting at 1. (I don’t remember.) AFAICS you need it to start at 1 and ramp down to -1 in the middle of the cycle. I’m thinking this because the A chain will jump at the beginning of the cycle, so you want A to be silent at that moment.

CombC will “ring” a little after the delay time jumps, but the LFTri will be at the silent point for only an instant. Maybe that’s not enough.

hjh

ouch, even after making some plots with weird outcome it was not coming to my mind. The Helpfile says for an initial phase for the LFSaw of 0 radians you need to set phase to 1; I think the values in the plot should be correct. LFTri starts at 1.

(
{
	var saw0 = LFSaw.ar(440, 1); // initial phase of 0 radians
	var saw90 = LFSaw.ar(440, 0); // initial phase of pi radians
	var tri90 = LFTri.ar(440, 1); // initial phase of pi radians

	[saw0, saw90, tri90];
}.plot;
)

Now there are no hard resets in the outcome but it still sounds kinds of rough at the crossing point. probably because of the flanger design with CombC and the ringing time.

(
SynthDef(\flangerSync, {

	var rate = \rate.kr(0.1);
	var minDelay = \minDel.kr(0.004);
	var maxDelay = \maxDel.kr(0.01);
	var decayTime = \decayTime.kr(1);
	var sig, sigA, sigB, lfoA, lfoB, combA, combB, delayTime;

	lfoA = LFSaw.ar(rate, 1, (maxDelay * 0.5) - minDelay, (maxDelay * 0.5) + minDelay); // initial phase of 0 radians
	lfoB = LFSaw.ar(rate, 0, (maxDelay * 0.5) - minDelay, (maxDelay * 0.5) + minDelay); // initial phase of pi radians

	sig = Saw.ar(100);

	combA = CombC.ar(sig, maxDelay, lfoA, decayTime);
	combB = CombC.ar(sig, maxDelay, lfoB, decayTime);

	sig = XFade2.ar(combA, combB, LFTri.kr(rate, 1)); // initial phase of pi radians

	sig = Pan2.ar(sig, \pan.kr(0), \amp.kr(0.1));
	//sig = SafetyLimiter.ar(sig);
	Out.ar(\out.kr(0), sig);
}).play;
)