BLIT headache

[juha_p]

...
Dirichlet kernels are "periodic sincs"

OK, I see, but..., when I ran couple comparison examples from Matlab site, I got results showing the equality is not fully perfect.

QUOTE: The Dirichlet and sinc functions are related by DN(πx)=sinc(Nx/2)/sinc(x/2). Relationship for N=6 while avoiding indeterminate expressions by specifying that the ratio of sinc functions is (−1)k(N−1) for x=2k, where k is an integer.

Code: Select all

% needed for Octave ----------
pkg load signal
% -----------------------------------
xmax = 4;
x = linspace(-xmax,xmax,1001)';
N = 6;
yd = diric(x*pi,N);
ys = sinc(N*x/2)./sinc(x/2);
ys(~mod(x,2)) = (-1).^(x(~mod(x,2))/2*(N-1));
figure(1);
plot(x,yd,ys)
legend('D_6(x*pi)', 'sinc(6*x/2) / sinc(x/2)', 'location', 'southwest');
title('The Dirichlet and sinc functions are related by DN(πx)=sinc(Nx/2)/sinc(x/2)... ')

N = 13;
yd = diric(x*pi,N);
ys = sinc(N*x/2)./sinc(x/2);
ys(~mod(x,2)) = (-1).^(x(~mod(x,2))/2*(N-1));
figure(2);
plot(x,yd, ys)
title('when N=13.')
legend('D_{13}(x*pi)', 'sinc(13*x/2) / sinc(x/2)', 'location', 'southwest')

dirsinc1.png

Zoomed:

dirsinc1b.png

dirsinc2.png

Zoomed:

dirsinc2b.png

[Z1202]

...
Dirichlet kernels are "periodic sincs"

OK, I see, but..., when I ran couple comparison examples from Matlab site, I got results showing the equality is not fully perfect.

I'm not sure what exactly you compared, but the equations listed there are not the ones describing the "periodic sinc" feature. If anything, then that one would be something like (very roughly, the functions still need to be properly scaled!)
D_n(t) = sum_i sinc(t+i)
That is Dirichlet kernel is a sum of unwindowed sincs (and hence a perfect BLIT)
But I was simply referring to the fact that if you take a completely flat spectrum, it's Dirac delta in continuous spectra and impulse train in discrete spectra. By bandlimiting those you respectively obtain sinc and D_n.

[elena2]

Thanks for the replies guys. This is really getting a more complex subject than I could imagine initially... anyway after much brainstorming and experimenting, I found a solution which apparently gives me a perfect aliasing free spectrum. It is about using DSF (discrete summation formulas).
1/n Sigma[k=1...n] cos(2PIkx) =
sin((2n+1)PIkx)/sin(PIkx)/2n
by using an integer n such that the highest harmonic is at most and not over Nyquist. It works wonders. To increase performance I have implemented it as a tabulated form. Since ultimately it is about the ratio of an integer multiple of a fundamental frequency and a fundamental frequency, I used a huge (262144 seems ok) sine table, large enough not to need interpolation (to avoid too many (int) casts to extract the fractional part) The denominator component is computed by addressing the table from a cyclical normalized running phase (i.e a plain sin oscillator), and the numerator component is computed by multiplying the table index by 2n+1 wrapping to table size by a simple And operation since the table size is a power of 2. I am quite happy with the result by now. There process is stable and can be freely modulated without producing artifacts and the spectrum is always a perfect comb.

Btw my DWG actually uses a crossfading between two taps when the pitch is modulated, but that is not enough to warrant stability to the one pole IIR all pass interpolator on the long run and with a feedback of unity to be used as a BLIT oscillator.

Btw I think I have understood more intuitively the flaw with my initial approach of placing every impulse shifted by a proper amount by a fractional delay filter: the impulse gets no leading ringing, only trailing of course, so the result is not a sinc but just half of it...

[Z1202]

I found a solution which apparently gives me a perfect aliasing free spectrum. It is about using DSF (discrete summation formulas).
1/n Sigma[k=1...n] cos(2PIkx) =
sin((2n+1)PIkx)/sin(PIkx)/2n

That's the Dirichlet kernel (possibly scaled)

the impulse gets no leading ringing, only trailing of course, so the result is not a sinc but just half of it...

That per se is not necessarily a problem. Depending on what kind of allpass you are using. But it may be.

[elena2]

I didn't know that was a Dirichlet kernel... I knew of it in a different context... well ok good to know !

[2DaT]

Thanks for the replies guys. This is really getting a more complex subject than I could imagine initially... anyway after much brainstorming and experimenting, I found a solution which apparently gives me a perfect aliasing free spectrum. It is about using DSF (discrete summation formulas).
1/n Sigma[k=1...n] cos(2PIkx) =
sin((2n+1)PIkx)/sin(PIkx)/2n
by using an integer n such that the highest harmonic is at most and not over Nyquist. It works wonders. To increase performance I have implemented it as a tabulated form. Since ultimately it is about the ratio of an integer multiple of a fundamental frequency and a fundamental frequency, I used a huge (262144 seems ok) sine table, large enough not to need interpolation (to avoid too many (int) casts to extract the fractional part) The denominator component is computed by addressing the table from a cyclical normalized running phase (i.e a plain sin oscillator), and the numerator component is computed by multiplying the table index by 2n+1 wrapping to table size by a simple And operation since the table size is a power of 2. I am quite happy with the result by now. There process is stable and can be freely modulated without producing artifacts and the spectrum is always a perfect comb.

DSF can cause clicks or other artifacts when number of harmonics change.

I'd rather use an approach similar to BLEPs - insert an "impulse" per phase wrap. It has a downside that it costs more the higher the frequency, but if it's implemented correctly it can be quite fast.

[elena2]

I have still to detect clicks and artifacts, by now it *seems* there are not any audible ones but I am still testing. The issue with BLEP from what I have understood is that it is an approximation, so the spectrum won't be a perfect comb as with the DSF method...

[Z1202]

I have still to detect clicks and artifacts, by now it *seems* there are not any audible ones but I am still testing. The issue with BLEP from what I have understood is that it is an approximation, so the spectrum won't be a perfect comb as with the DSF method...

OTOH the issue with BLIT is that it has problems if the frequency is not constant (and clicks from switching between different order kernels is probably the least of the problems).

[mystran]

I have still to detect clicks and artifacts, by now it *seems* there are not any audible ones but I am still testing. The issue with BLEP from what I have understood is that it is an approximation, so the spectrum won't be a perfect comb as with the DSF method...

Yes, they are approximations, but the approximation is based on standard interpolation theory (same stuff your audio interface uses to actual send the signal to your speakers) and you can make your approximations as good as you deem necessary. Eventually you'll hit floating point noise-floor anyway.

Yes, they need more CPU as the fundamental gets higher, but realistically with a decent implementation that usually only starts to become a significant issue once your fundamental starts getting into the megahertz range. In more typical use, they are usually cheaper than calling sin() twice per sample.