obscure antialiased saw code [Ancient thread bumped; see page 3]

DSP, Plug-in and Host development discussion.
KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Mon Feb 20, 2006 3:17 am

aciddose wrote:i have been unable to find the function to calculate the correct gain factor based upon frequency. the function most likely is more complex than the generator itself.
Here is the code as provided by thierry rochebois, I translated and added a few comments.

Code: Select all

/* clair.c         Examen Partiel 2b
   T.Rochebois
   02/03/98
*/
#include <stdio.h>
#include <math.h>
main()
{
  double phase=0,dphase,freq,compensation;
  double aw0=0,aw1=0,ax0=0,ax1=0,ay0=0,ay1=0,az0=0,az1=0,sortie;
  short aout;
  int sr=44100;       //sample rate (Hz)
  double f_debut=55.0;//start freq (Hz)
  double f_fin=sr/6.0;//end freq (Hz)
  double octaves=log(f_fin/f_debut)/log(2.0);
  double duree=50.0;  //duration (s)
  int i;
  FILE* f;
  f=fopen("saw.pcm","wb");
  for(i=0;i<duree*sr;i++)
  {
    //exponential frequency sweep
    //Can be replaced by anything you like.
    freq=f_debut*pow(2.0,octaves*i/(duree*sr));
    dphase=freq*(2.0/sr);     //normalised phase increment
    phase+=dphase;            //phase incrementation
    if(phase>1.0) phase-=2.0; //phase wrapping (-1,+1)

    //polynomial calculation (extended continuity at -1 +1)
    //        7       1  3    1   5
    //P(x) = --- x - -- x  + --- x
    //       360     36      120
    aw0=phase*(7.0/360.0 + phase*phase*(-1/36.0 + phase*phase*(1/120.0)));
    // quad differentiation (first order high pass filters)
    ax0=aw1-aw0; ay0=ax1-ax0; az0=ay1-ay0; sortie=az1-az0;
    //compensation of the attenuation of the quad differentiator
    //this can be calculated at "control rate" and linearly
    //interpolated at sample rate.
    compensation=1.0/(dphase*dphase*dphase*dphase);
    // compensation and output
    aout=(short)(15000.0*compensation*sortie);
    fwrite(&aout,1,2,f);
    //old memories of differentiators
    aw1=aw0; ax1=ax0; ay1=ay0; az1=az0;
  }
  fclose(f);
}
The compensation function is quite simple.
Last edited by Paul Sernine on Mon Jun 19, 2006 3:26 am, edited 1 time in total.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Mon Feb 20, 2006 3:27 am

Fire Sledge - Ohm Force wrote:I played a bit with this oscillator and it has a big drawback. It may click when pitch is changed. So it is required to recompute the whole oscillator state under some conditions.
-- Laurent
That's right. I believe it is a drawback of high pass filters.
I used this algo for PWM (done with two dephased sawteeth), the fast phase modulation (modulation by a fast Attack of an ADSR) induced a fast freq modulation and a click. As it was a transitory artefact (and not periodic) I simply clipped the output of the oscillator.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

User avatar
KVRAF
12422 posts since 7 Dec, 2004

Post Mon Feb 20, 2006 10:18 am

i still dont see how 24046326 relates to 1/(cps^4).

the compensation required for the differentiation is simple, i was just having problems finding the relation between the numbers actually used in the original code and the math that is to be expected.

anyway it isnt important to me, i use this instead
http://xhip.cjb.net/temp/saw.wav

KVRist
121 posts since 2 Nov, 2000 from 404 - Not found

Post Mon Feb 20, 2006 10:52 am

For the compensation, I found the same formula: gain = 1 / (360 * step^4), the 4th derivate of the polynomial f(x) = 3*x^5 - 10*x^3 + 7 being f""(x) = 360*x.

I also found an empirical compensation for the volume loss in high freq, but valid up to nyquist/2 : gain /= 1 - step * 1.15;

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Tue Feb 21, 2006 1:04 am

Just an idea. Maybe if we do the compensation before the differentiations part of the clicks can be removed. I will give it a try asap (i am very busy doing my boring job).
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist
121 posts since 2 Nov, 2000 from 404 - Not found

Post Tue Feb 21, 2006 3:27 am

Paul Sernine wrote:Maybe if we do the compensation before the differentiations part of the clicks can be removed.
I doubt it will do something, because clicks come mainly from the step length change. Thus, states of the derivates don't correspound to the value they should have after the frequency change. The difference is very small, but is extremely amplified with the successive differenciations.

-- Laurent

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Thu Feb 23, 2006 12:40 am

Fire Sledge - Ohm Force wrote:
Paul Sernine wrote:Maybe if we do the compensation before the differentiations part of the clicks can be removed.
I doubt it will do something, because clicks come mainly from the step length change. Thus, states of the derivates don't correspound to the value they should have after the frequency change. The difference is very small, but is extremely amplified with the successive differenciations.

-- Laurent
I tried very empirically to improve the results by changing the way the compensation is done. Here is the way I obtained the best results:

Code: Select all

/* clair.c         Examen Partiel 2b
   T.Rochebois
   02/03/98
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
main()
{
  double phase=0,dphase,freq,comp;
  double aw0=0,aw1=0,ax0=0,ax1=0,ay0=0,ay1=0,az0=0,az1=0,fout;
  short aout;
  int sr=44100;
  double f_min=55.0;
  double f_max=sr/6.0;
  double octaves=log(f_max/f_min)/log(2.0);
  double duree=50.0;
  double freqf=55.0,freqf2=55.0,freqf3=55.0,freqf4=55.0;

  int i;
  FILE* f;
  f=fopen("saw.pcm","wb");
  for(i=0;i<duree*sr;i++)
  {
    if((i%3000)==0)
      freq=f_min*pow(2.0,octaves*(rand()*(1.0/RAND_MAX)));

    freqf +=0.05*(freq  -freqf);
    freqf2+=0.05*(freqf -freqf2);
    freqf3+=0.05*(freqf2-freqf3);
    freqf4+=0.05*(freqf3-freqf4);

    dphase=freqf4*(2.0/sr);
    phase+=dphase;
    if(phase>=1.0) phase-=2.0;
    comp=1.0/dphase;
    aw0=phase*(7.0/360.0 + phase*phase*(-1/36.0 + phase*phase*(1/120.0)));
    ax0=comp*(aw1-aw0); ay0=comp*(ax1-ax0); az0=comp*(ay1-ay0); fout=comp*(az1-az0);
    fout=10000.0*fout;
    if(fout>32000) fout=32000;
    if(fout<-32000) fout=-32000;
    aout=(short)(fout);
    fwrite(&aout,1,2,f);
    aw1=aw0; ax1=ax0; ay1=ay0; az1=az0;
 }
  fclose(f);
}
It is still not perfect, just a little better. But I believe that there is a way to do it right.
I will try to be less empirical (and less lazy) and go analytic to find the correct way.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Mon Feb 27, 2006 2:04 am

I played with very fast and deep modulations this week end, without clicks :) .

Code: Select all

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
main()
{
  double phase=0,dphase,freq,comp;
  double aw0=0,aw1=0,ax0=0,ax1=0,ay0=0,ay1=0,az0=0,az1=0,fout;
  short aout;
  int sr=44100;
  double duree=10.0,f_min=100.0,f_max=2000.0; // sweep
  double lfospeed=16, lfodepth=6; //mod:16cps on 6=[-3,+3] octaves
  double octaves=log(f_max/f_min)/log(2.0);
  int i;
  FILE* f;
  f=fopen("saw.pcm","wb");
  for(i=0;i<duree*sr;i++)
  {
    freq=f_min*pow(2.0,octaves*i/(duree*sr)+0.5*lfodepth*sin(i*6.28*lfospeed/sr));
    dphase=freq*(2.0/sr);
    phase+=dphase;
    if(phase>=1.0) phase-=2.0;
    aw0=phase*(7.0/360.0 + phase*phase*(-1/36.0 + phase*phase*(1/120.0)));
    comp=1/dphase;
    ax0=(aw1-aw0)*comp;
    ay0=(ax1-ax0)*comp;
    az0=(ay1-ay0)*comp;
    fout=(az1-az0)*comp;

    fout=10000.0*fout;
    if(fout>32000) fout=32000;
    if(fout<-32000) fout=-32000;
    aout=(short)(fout);
    fwrite(&aout,1,2,f);
    aw1=aw0; ax1=ax0; ay1=ay0; az1=az0;
  }
  fclose(f);
}
I think that clicks appear when the derivatives of the frequency are not continuous enough.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRian
646 posts since 18 Feb, 2006 from California

Post Mon Feb 27, 2006 12:47 pm

in the 2k+ range of this algorithm (paul's version) there is an artifact that sounds like some kind of fm. What causes this?

here's a 2k-5k sweep, parodon slow server
http://music.calarts.edu/~asomers/saw_a ... s.pcm.aiff

would this be a very bad real-time saw generator?
aciddose wrote:you empty headed animal food trough wiper, your mother was a hamster and your father smelt of elderberries!

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Wed Mar 01, 2006 1:37 am

asomers wrote:in the 2k+ range of this algorithm (paul's version) there is an artifact that sounds like some kind of fm. What causes this?

here's a 2k-5k sweep, pardon slow server
http://music.calarts.edu/~asomers/saw_a ... s.pcm.aiff
Your aiff file is a stereo file ! The input file is a raw mono file, I think that's why it is aliased (hence the FM character :wink: ). Just import it as a mono file or double the fwrite and import it as a stereo file :

Code: Select all

...
aout=(short)(fout); 
fwrite(&aout,1,2,f);
fwrite(&aout,1,2,f);
...
asomers wrote: would this be a very bad real-time saw generator?
I don't think so. I think that it can be very usefull even if it is rather trivial to say and repeat and repeat and repeat that minBLET (or the "corrective grains" variant) is better.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist
121 posts since 2 Nov, 2000 from 404 - Not found

Post Wed Mar 01, 2006 1:49 am

Indeed it works well ! I fixed a bit my version, i did a post-derivation incrementation of the phase instead of pre- and it gave less tolerance to frequency change. Now it works better.

KVRist
101 posts since 25 Aug, 2004

Post Wed Mar 01, 2006 3:02 am

Off topic :

Paul Sernine??? Seriously?? Or are you a Leblanc fan?

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Wed Mar 01, 2006 4:42 am

gonzo wrote:Off topic :

Paul Sernine??? Seriously?? Or are you a Leblanc fan?
Sure, I am (a Leblanc fan) :D
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Wed Mar 01, 2006 4:51 am

Fire Sledge - Ohm Force wrote:Indeed it works well ! I fixed a bit my version, i did a post-derivation incrementation of the phase instead of pre- and it gave less tolerance to frequency change. Now it works better.
Just for fun, I improved it so that it is even more stable. The derivators are sort-of time delayed by half a sample at each stage.

Code: Select all

//Perfectionniste
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
main()
{
  double phase=0,dphase,freq;
  double aw0=0,aw1=0,ax0=0,ax1=0,ay0=0,ay1=0,az0=0,az1=0,fout;
  short aout;
  int sr=44100;
  double duree=50.0,f_min=50.0,f_max=3000.0; // sweep
  double lfospeed=40, lfodepth=8; // modulation 40 cps , 8 octaves
  double octaves=log(f_max/f_min)/log(2.0);
  double dpw0=1,dpw1=1,dpx0=1,dpx1=1,dpy0=1,dpy1=1,dpz0=1,dpz1=1;
  int i;
  FILE* f;
  f=fopen("saw.pcm","wb");
  for(i=0;i<duree*sr;i++)
  {
    freq=f_min*pow(2.0,octaves*i/(duree*sr)+0.5*lfodepth*sin(i*6.28*lfospeed/sr));

    dphase=freq*(2.0/sr);
    phase+=dphase;
    if(phase>=1.0) phase-=2.0;
    aw0=phase*(7.0/360.0 + phase*phase*(-1/36.0 + phase*phase*(1/120.0)));

    //ax0 is calculated on the short integration segment dpw0=dphase.
    dpw0=dphase;
    ax0=(aw1-aw0)/dpw0;

    //ay0 is calculated on the short integration segment dpx0.
    dpx0=0.5*(dpw0+dpw1);
    ay0=(ax1-ax0)/dpx0;

    //az0 is calculated on the short integration segment dpy0.
    dpy0=0.5*(dpx0+dpx1);
    az0=(ay1-ay0)/dpy0;

    //fout is calculated on the short integration segment dpz0.
    dpz0=0.5*(dpy0+dpy1);
    fout=(az1-az0)/dpz0;

    fout=10000.0*fout;
    if(fout>32000) fout=32000;
    if(fout<-32000) fout=-32000;
    aout=(short)(fout);
    fwrite(&aout,1,2,f);
    aw1=aw0;    ax1=ax0;   ay1=ay0;   az1=az0;
    dpw1=dpw0; dpx1=dpx0; dpy1=dpy0; dpz1=dpz0;
 }
  fclose(f);
}
Last edited by Paul Sernine on Wed Mar 01, 2006 9:04 am, edited 2 times in total.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

KVRist

Topic Starter

71 posts since 14 Feb, 2006 from Variable

Post Wed Mar 01, 2006 5:44 am

Comments : Time domain interpretation of Thierry Rochebois' polynomial algo.
The four differenciators (as opposed to derivators) are interpreted as four "means" on time segments dpw0,dpx0,dpy0,dpz0. Each "mean" is an integral divided by the width of the integration segment (the delta phases). This way, I guess that the differentiators can be interpreted as a series of very-small-box filters.
Last edited by Paul Sernine on Mon Jun 19, 2006 3:28 am, edited 1 time in total.
Science is the belief in the ignorance of the experts.
Richard P. F
EYNMAN

Return to “DSP and Plug-in Development”