Hi, I have already asked on dsp.stackexchange.com without answers, can you help me ?
In a Kraplus Strong algorithm we experience detuning because of the LP filter in the feedback line. I want to compensate for this detuning but to do that I need to know the formula relating filtering factor k and resulting pitch drift. I failed to find a suitable formula, the underlying math gets too complex for me. In short, if we call n the (non integer) delay lenght corrisponding to the required pitch, and k the filter factor (lp filter in form of y=x+k(y[-1]-x) ), I need to know the precise math formula linking k and n' where n' is >n and corresponds to the detuning I get, so I can use a compensated n to get rid of detuning in function of k. The only things I discovered of this mystery function is that it passes by 0 for k=0, by 1 for k=0.5 and likely goes to infinity for k=1. But it is not k/(1-k) even if it is a fair approximation, and it is not even (k/(1-k))^n or ^1/n either.
Can you help me please ? Thank in advance.
How to compensate for LP detuning in Karplus Strong
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
Right, so the phase-shift that results in the detuning is frequency dependent (unless you use a linear-phase damping filter), so each harmonic is actually detuned by a different amount. On the bright side (since fixing them all is a lot more complicated), as long as the fundamental is correctly tuned, slightly inharmonic partials are usually perceived mostly as an alternation of the tone and few real-life instruments are perfectly harmonic either.
To compensate for the tuning of the fundamental, you need to compute the phase-shift at the fundamental frequency. Unfortunately this requires at least a bit of complex arithmetic. For a one-pole low-pass such as your's, you have the transfer function (1-k)/(1-k*z^-1) [edit: fixed my original mistake with signs, sorry about that], where we substitute z^-1=exp(-i*2*pi*f/fs) where "f" is the frequency of interest (eg. the desired fundamental) and "fs" is the sampling rate (and "i" is the imaginary unit), then compute the complex argument of the result (which gives us the phase-shift; since the filter is first order, there is no need to worry about wrap around). Then to compute the actual time-delay in samples, divide the computed phase by 2*pi*f/fs.
Unfortunately there really isn't any short-cut to this math.
Strictly speaking, unless you are using a decent-quality sinc-interpolator to implement your fractional delays, your interpolation probably adds a bit of detune as well, especially at very high frequencies. This is typically much less of a problem than the detuning from a damping filter (and requires somewhat more complicated math to work around), but I still wanted to mention it, in case you find that even after compensating for the damping, your fundamental is slightly off at very high frequencies.
To compensate for the tuning of the fundamental, you need to compute the phase-shift at the fundamental frequency. Unfortunately this requires at least a bit of complex arithmetic. For a one-pole low-pass such as your's, you have the transfer function (1-k)/(1-k*z^-1) [edit: fixed my original mistake with signs, sorry about that], where we substitute z^-1=exp(-i*2*pi*f/fs) where "f" is the frequency of interest (eg. the desired fundamental) and "fs" is the sampling rate (and "i" is the imaginary unit), then compute the complex argument of the result (which gives us the phase-shift; since the filter is first order, there is no need to worry about wrap around). Then to compute the actual time-delay in samples, divide the computed phase by 2*pi*f/fs.
Unfortunately there really isn't any short-cut to this math.
Strictly speaking, unless you are using a decent-quality sinc-interpolator to implement your fractional delays, your interpolation probably adds a bit of detune as well, especially at very high frequencies. This is typically much less of a problem than the detuning from a damping filter (and requires somewhat more complicated math to work around), but I still wanted to mention it, in case you find that even after compensating for the damping, your fundamental is slightly off at very high frequencies.
Last edited by mystran on Thu Feb 13, 2020 2:49 am, edited 1 time in total.
-
- KVRian
- 1000 posts since 1 Dec, 2004
You need to compute either the "group delay" of your filter, or the phase offset caused by your filter at the desired root pitch.
Unfortunately, I can't find any good info about calculating group delay which aren't basically walls of math.
You will also need to be able to have a delay length that isn't integer (otherwise high notes will never be in tune). See this classic paper for a way that works well with physical modelling (1-pole allpass filter): https://ccrma.stanford.edu/~dattorro/Ef ... nPart2.pdf
Unfortunately, I can't find any good info about calculating group delay which aren't basically walls of math.
You will also need to be able to have a delay length that isn't integer (otherwise high notes will never be in tune). See this classic paper for a way that works well with physical modelling (1-pole allpass filter): https://ccrma.stanford.edu/~dattorro/Ef ... nPart2.pdf
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Thanks a lot for your precise answer. I will try to work out your math. To me it still remains a mystery why my formula k/(1-k), by accurate numerical analysis carried with fft, gives a so precise estimation of the delta n only up to k=0.9 or so and then diverges quite abruptly. But it is not in function of the fundamental frequency so clearly it can not be precise.
As for interpolation I am using simple linear interpolation reading the circular buffer position taking the fractional part as sub sample index for linear itp. This approach always worked great for me with no audible tuning defects, the only (tolerable) drawback being that with fundamentals corresponding to zero or close to zero fractional part of n, the sound is crispier (less hf filtering) but that is obvious. Otoh, implementing more complex interpolations like a sinc or all pass filter as suggested by some authors seems an overkill to me.
As for interpolation I am using simple linear interpolation reading the circular buffer position taking the fractional part as sub sample index for linear itp. This approach always worked great for me with no audible tuning defects, the only (tolerable) drawback being that with fundamentals corresponding to zero or close to zero fractional part of n, the sound is crispier (less hf filtering) but that is obvious. Otoh, implementing more complex interpolations like a sinc or all pass filter as suggested by some authors seems an overkill to me.
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
For resonant frequencies, you need phase delay (which is simply phase-shift times wavelength) rather than group delay (which would be the derivative of phase response over frequencies). The latter is a measure of the delay on the envelope of the signal, but it's irrelevant for computing the resonant frequencies themselves, which happen precisely at those frequencies where the total phase-shift (including the delay line) is a multiple of 2*pi.
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Mystran, it was an effort but at the end I could come to the following formula for delta n (i.e the drift of n in function of k) basing on your precious suggestion:
Delta n=atan(2Pi*k/(n*(1-k)))*n/2Pi
...which is in fact VERY CLOSE to my k/(1-k) estimation for k<0.9.
It is enough to subtract this value from.the original n (delay lenght) at it apparently does its job
Delta n=atan(2Pi*k/(n*(1-k)))*n/2Pi
...which is in fact VERY CLOSE to my k/(1-k) estimation for k<0.9.
It is enough to subtract this value from.the original n (delay lenght) at it apparently does its job
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
First of all, looks like I made a small mistake with the transfer function, it should obviously be (1-k)/(1-k*z^-1). I'm terribly sorry if you wasted time due to this, sometimes you just get these small errors when you type things quickly.
atan2(0,1-k)-
atan2((sin((2*%pi*f)/fs)*k)/sqrt((1-cos((2*%pi*f)/fs)*k)^2+sin((2*%pi*f)/fs)^2*k^2)
,(1-cos((2*%pi*f)/fs)*k)/sqrt((1-cos((2*%pi*f)/fs)*k)^2+sin((2*%pi*f)/fs)^2*k^2))
This looks a bit scarier than it actually is, because we know the second argument to both atan2() functions is always non-zero; the first atan2() is then always zero and we can write the second one using regular atan() with division, so the square roots cancel out and we get:
atan(sin(k*2*pi*f/fs)/(1-k*cos(2*pi*f/fs))
= atan(sin(k*2*pi/n) / (1-k*cos(2*pi/n)), assuming n=fs/f is the desired total delay in samples.
This is the exact result for the phase (well, assuming I didn't make any more mistakes).
However, for low-enough frequencies, we could make some further approximations (although these will be quite poor once you start getting into the kHz range with typical sampling rates). For small enough x, sin(x) and atan(x) are both approximately equal to x and cos(x) is approximately equal to 1.
If we take all of these short-cuts, we get (k*2*pi/n)/(1-k). Again, this is the phase, so to get the delay adjustment, multiply by n/(2*pi). The formula you arrived at then is what you get if you make the sin() and cos() approximations and leave the atan as-is. That said, I would probably just compute the sines and cosines properly (or at least use slightly better approximations).
Anyway, let me try to explain this formula of yours. From the outline I gave you (with the stupid transfer function mistake fixed, lol), Maxima (the computer algebra system I usually use, the percentages are due to copy-paste from Maxima) gives:
atan2(0,1-k)-
atan2((sin((2*%pi*f)/fs)*k)/sqrt((1-cos((2*%pi*f)/fs)*k)^2+sin((2*%pi*f)/fs)^2*k^2)
,(1-cos((2*%pi*f)/fs)*k)/sqrt((1-cos((2*%pi*f)/fs)*k)^2+sin((2*%pi*f)/fs)^2*k^2))
This looks a bit scarier than it actually is, because we know the second argument to both atan2() functions is always non-zero; the first atan2() is then always zero and we can write the second one using regular atan() with division, so the square roots cancel out and we get:
atan(sin(k*2*pi*f/fs)/(1-k*cos(2*pi*f/fs))
= atan(sin(k*2*pi/n) / (1-k*cos(2*pi/n)), assuming n=fs/f is the desired total delay in samples.
This is the exact result for the phase (well, assuming I didn't make any more mistakes).
However, for low-enough frequencies, we could make some further approximations (although these will be quite poor once you start getting into the kHz range with typical sampling rates). For small enough x, sin(x) and atan(x) are both approximately equal to x and cos(x) is approximately equal to 1.
If we take all of these short-cuts, we get (k*2*pi/n)/(1-k). Again, this is the phase, so to get the delay adjustment, multiply by n/(2*pi). The formula you arrived at then is what you get if you make the sin() and cos() approximations and leave the atan as-is. That said, I would probably just compute the sines and cosines properly (or at least use slightly better approximations).
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Yeah is is a mess, but indeed I found the approximation does its job nicely. I didn' t even notice the typo don' t mind since I assumed the correct transfer formula implicitely. Since k and n are allowed to change on a sample rate basis, I would be better using tables though for actual computation.. I will now test well both versions though to decide if it is (acoustically) worth the trouble. I agree that to simulate natural sounds like plucked strings one should allow some imprecision for better realism. But the detuning caused by damping was just too big and intolerable. Now at least I have a formula to fix it. Thanks again for your patience
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
I would guess (though you really need to measure these things) than at least on a desktop CPU, computing it on the fly is probably going to be faster than reading it from a lookup table.
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Maybe you' re even right since in general accessing memory is slower than doing inline computations but it depends (on my i7 I am using sin cos tables because after testing they still result faster than having to call the math lib functions.) But really it depends...
Anyway I have now tried enough and I can say that my approximation (without the sin and cos inside the atan) is pretty perfect. The error is neither audible not even easily measureable so I am happy
Since we touched the interpolation topic, I must specify a thing here in case somebody else may find it useful in future.
Linear interpolation for fractional delays is unbeatable for speed and simplicity but has the known drawback of filtering high frequencies depending on fr (fr=delay_time_real - (int)delay_time_real). When fr=0 there is no filtering, while at 0.5 there is max filtering and Nyquist is gone - for obvious reasons.
Many authors suggest using better solutions like an all-pass filter or even sinc interpolation, which are still expensive and complex.
I invented an evil trick which after deep testing seems to work perfectly and likely offers an almost identical result of actual sinc interpolation. A really genial invention allow me to say, so I wanted to share it.
After plain linear interpolation i.e sample=buffer[n]*(1-fr) + buffer[n+1]*fr,
a simple comb filter centered at Nyquist (i.e with delay=1 sample) and negative feedback factor b=-fr' ^2 is applied, where fr' is computed as follows:
fr'=2*fr;
if(fr' >1) fr'=2-fr';
so a comb filter like this:
y[t]=x[t] -fr'^2 * (y[t-1]-x[t])
which has the effect of compensating perfectly the hf damping centered at Nyquist of the linear interpolation in function of fr.
Analysis of such linear interpolation plus compensating comb filter shows flat spectrum response at the analyzer upto the very Nyquist frequency for all values of fr. The effect on a square wave is the tiny ringing at edges one would expect from sinc interpolation.
I suspect I discovered a very smart math shortcut to avoid all the computations usually needed for sinc interpolation, other than having elegantly solved a very old and known issue with fractional delays...
Anyway I have now tried enough and I can say that my approximation (without the sin and cos inside the atan) is pretty perfect. The error is neither audible not even easily measureable so I am happy
Since we touched the interpolation topic, I must specify a thing here in case somebody else may find it useful in future.
Linear interpolation for fractional delays is unbeatable for speed and simplicity but has the known drawback of filtering high frequencies depending on fr (fr=delay_time_real - (int)delay_time_real). When fr=0 there is no filtering, while at 0.5 there is max filtering and Nyquist is gone - for obvious reasons.
Many authors suggest using better solutions like an all-pass filter or even sinc interpolation, which are still expensive and complex.
I invented an evil trick which after deep testing seems to work perfectly and likely offers an almost identical result of actual sinc interpolation. A really genial invention allow me to say, so I wanted to share it.
After plain linear interpolation i.e sample=buffer[n]*(1-fr) + buffer[n+1]*fr,
a simple comb filter centered at Nyquist (i.e with delay=1 sample) and negative feedback factor b=-fr' ^2 is applied, where fr' is computed as follows:
fr'=2*fr;
if(fr' >1) fr'=2-fr';
so a comb filter like this:
y[t]=x[t] -fr'^2 * (y[t-1]-x[t])
which has the effect of compensating perfectly the hf damping centered at Nyquist of the linear interpolation in function of fr.
Analysis of such linear interpolation plus compensating comb filter shows flat spectrum response at the analyzer upto the very Nyquist frequency for all values of fr. The effect on a square wave is the tiny ringing at edges one would expect from sinc interpolation.
I suspect I discovered a very smart math shortcut to avoid all the computations usually needed for sinc interpolation, other than having elegantly solved a very old and known issue with fractional delays...
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
The standard math functions are typically somewhat slow, because they have to handle all the possible special cases, so using direct approximations (that ignore corner cases that are irrelevant for audio) is usually faster. Additionally, you can typically get large speedups by processing using SIMD.elena2 wrote: ↑Thu Feb 13, 2020 8:28 pm Maybe you' re even right since in general accessing memory is slower than doing inline computations but it depends (on my i7 I am using sin cos tables because after testing they still result faster than having to call the math lib functions.) But really it depends...
While sinc-interpolation is more expensive, it also has a problem with feedback loops in the sense that the (unavoidable) latency puts a limit on the minimum delay length (eg. if you have 16 samples of latency, then you can't have a loop shorter than 16 samples, because that would require negative delay). One can reduce this problem by oversampling, but then it gets even more expensive. I would usually suggest "cubic-interpolation" (ie. Catmull-Rom hermite splines; Wikipedia has the math you need) as the next step up from linear. This requires 4 points (rather than 2, the "latency" is only 1 extra sample though, so not really a problem in most cases) and slightly more computation, but the results are also much better.Linear interpolation for fractional delays is unbeatable for speed and simplicity but has the known drawback of filtering high frequencies depending on fr (fr=delay_time_real - (int)delay_time_real). When fr=0 there is no filtering, while at 0.5 there is max filtering and Nyquist is gone - for obvious reasons.
Many authors suggest using better solutions like an all-pass filter or even sinc interpolation, which are still expensive and complex.
Filters with single-sample delays are not normally called "comb filters" (ie. the "comb" name comes from the shape of the response you get when you have a longer delay, since it has "comb-like" spikes at regular intervals) but rather just "filters" and this is really just an example of post-equalisation. This is a valid approach for (approximately) flattening the spectrum, but it will also boost the interpolation artefacts. Additionally using slightly better interpolation (eg. cubic) might actually be cheaper. There are many possible solutions though.After plain linear interpolation i.e sample=buffer[n]*(1-fr) + buffer[n+1]*fr,
a simple comb filter centered at Nyquist (i.e with delay=1 sample) and negative feedback factor b=-fr' ^2 is applied, where fr' is computed as follows:
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Yep I know it doesn' t look like a comb anymore but it is still basically one if you use a delay longer than one sample.
Anyway my trick is performing really well by now, pitch is reproduced with perfect continuity and I am not noticing or hearing any artifacts... if that eventually turns out the case I will try with cubic interpolation, I have my personal spline formulation with 4 points indeed which I have always been using for sample playing purposes
Anyway my trick is performing really well by now, pitch is reproduced with perfect continuity and I am not noticing or hearing any artifacts... if that eventually turns out the case I will try with cubic interpolation, I have my personal spline formulation with 4 points indeed which I have always been using for sample playing purposes
-
- KVRist
- 88 posts since 5 Jul, 2018 from Cambridge, UK
So, the concept of combining an FIR filter (the linear interpolation) with an "opposite" feedback filter is very similar to the construction of an all-pass filter.elena2 wrote: ↑Thu Feb 13, 2020 8:28 pm After plain linear interpolation i.e sample=buffer[n]*(1-fr) + buffer[n+1]*fr,
a simple comb filter centered at Nyquist (i.e with delay=1 sample) and negative feedback factor b=-fr' ^2 is applied, where fr' is computed as follows:
fr'=2*fr;
if(fr' >1) fr'=2-fr';
so a comb filter like this:
y[t]=x[t] -fr'^2 * (y[t-1]-x[t])
which has the effect of compensating perfectly the hf damping centered at Nyquist of the linear interpolation in function of fr.
This is definitely useful for inter-sample interpolation: https://ccrma.stanford.edu/~jos/pasp/Fi ... ation.htmlThere are a few issues I can see with your proposed equations: firstly, it looks like your filter is unstable at fr=0.5.
It also doesn't look (to me) like your feedback filter perfectly compensates for the linear-interpolation damping. When it comes down to it, I think you end up using fr' = (2*fr)^2 as an approximation of fr/(1 - fr) for the range 0 < fr < 0.5. So it's closer than not compensating at all, but I'd be surprised if it's exactly flat.
But also, the feedback filter is introducing its own delays! (Which will also vary by frequency, like the lowpass which was the original issue.) The construction in the first-order-allpass link above ends up splitting the delays between the two parts of the filter, and there are a few more subtleties from that as well.
When you checked your filter's behaviour, how did you verify it was "flat", and did you check the phase response (because it will no longer be linear, due to the feedback part)?
-
- KVRist
- Topic Starter
- 109 posts since 7 Jun, 2018
Well I am a coder first of all and not a mathematician lol... so I leave to somebody else the hyerogliphs to demonstrate what the combination of the two filters correspond to.. to me it looks like the result is pretty the same of a sinc interpolation, by wave analysis by naked eye. Often I achieve more of my goals by intuition that by filling papers of formulas... I have the linear itp process which creates a valley at Nyquist which cancels it completely for fr=0.5 as it is expected. So I compensate downstream by doing the opposite task: using a filter which "resonates" nyquist nothing at all for fr=0 and completely for fr=0.5. By some trials I discovered the quadratic relationship. I use my spectral analysis framework for analyzing the result. I fed a fractional delay coded as described with compensating filter with white noise and I analyzed the result with a 16384 fft frame size by varying fr between 0 and 1. The spectrum is FLAT, period. Also, I could detect no phase shift. I also fed the stuff with a square wave and with sine waves of varying frequency overimposed with a scope on the unfiltered signal. Ok for fr=0.5 the second filter is perhaps instable but for fr=0.5 nyquist is zero because the interpolator suppresses it.
People I am the author of many synths and of the spectral modules for Synthedit so I am not a stupid but one who definitely knows what she does. I have been coding audio and image processing algorithms since I was 14 or so... Since yesterday I am playing my synth with this trick implemented (pretty I was trying to perfection the tuned delay line of my Synlay synth) and everything sounds great. So... what else should I conclude ? That my trick works.
People I am the author of many synths and of the spectral modules for Synthedit so I am not a stupid but one who definitely knows what she does. I have been coding audio and image processing algorithms since I was 14 or so... Since yesterday I am playing my synth with this trick implemented (pretty I was trying to perfection the tuned delay line of my Synlay synth) and everything sounds great. So... what else should I conclude ? That my trick works.
-
- KVRAF
- 7890 posts since 12 Feb, 2006 from Helsinki, Finland
No need to get defensive, we are not trying to criticize (certainly not you as a person), but simply trying to help.elena2 wrote: ↑Fri Feb 14, 2020 5:40 pm Well I am a coder first of all and not a mathematician lol... so I leave to somebody else the hyerogliphs to demonstrate what the combination of the two filters correspond to.. [...] I have been coding audio and image processing algorithms since I was 14 or so... Since yesterday I am playing my synth with this trick implemented (pretty I was trying to perfection the tuned delay line of my Synlay synth) and everything sounds great. So... what else should I conclude ? That my trick works.
That said, don't be afraid of the math. It's scary at first, but it's worth it. At the end of the day, it isn't all that different from programming (and one could argue that when you use a compute algebra system to do most of the actual work for you, it's basically just another type of programming really), there is a bit of a learning curve at first (especially with DSP, where you kinda have to get used to working with complex numbers), but anyone capable of programming should be able to pickup the basics with just a little bit of effort.
