Beginner JS questions and help..

[saddle]

I used to program all the time, but I'm brand new to DSP stuff, and JS in general.

I immersed myself in FFT theory over the last few days, and I now have a headache! But I at least understand some of the principals. However in trying to decipher some of the math in the JS plugins that come with Reaper, I've come across some common math/functions that I am just curious about what they do.

For example from the Utility/Volume JS

@slider
adj1=2 ^ (slider1/6);
adj2=2 ^ (slider2/6);

What is the purpose of the 2^(slider1/6) ? It appears to convert the input (a db value) to a +- value, but what, why?

Are the samples a standard PCM value? Would you need to scale them before and after modification? I remember reading that the f(t) values had to be multiplied before use.

Just trying to figure out some of the basics so I can get a better understanding and try to code something.

The FFT functions to apply filters, eq, etc. I'm slowly finding and attempting to understand how to use, but any direction (and code snippits) on the general use of FFT functions to manipulate samples would be great!

I can't really find any tutorials on JS so far, and I've been looking through all 51 pages of this forum. But I'm sure I would know and understand the answer if I had run across it.

Thanks in advance for your help, patience and direction!
Saddle

[ashcat_lt]

The sample values are basically PCM values. They represent the fraction of the full output voltage that you want your soundcard to go to for that sample period. They usually run from -1 to 1, which represents 0dbfs.

The sliders in JSVolume are in db, not real values, or even real ratios. The first one sets the gain. If you want to add 6db, what you actually need to do is multiply the sample value by 2. Slider2 sets the limit level in dbfs. If you want it to limit at -6dbfs, what you need to do is limit the sample values to half what it would be at 0dbfs - so they will be between -0.5 and 0.5. Those statements sort of turn the db values into real factors or ratios, though they're not exactly accurate.

[SaulT]

ratio = 2^(decibel/6)

This is an approximation of the more proper

ratio = 10^(decibel/20)

It's a pretty decent approximation (within 1% accuracy) until you get to around -50 dB, and even then it's still pretty close (0.0031004 vs 0.0031623). In certain architectures, the first is a little faster. In JS, being interpreted, it doesn't really matter. Probably the only really important thing is to be consistent - either use one or the other, not both.

Anyway, if that needs a little further explanation, decibel is a ratio. Adding 1 dB is the same as multiplying by 1.122, subtracting 6 dB is about the same as dividing by 2. You can think of it as a shorthand, a way to think of audio in an easy linear fashion vs a more conceptually difficult logarithmic one. +1 dB is a small boost. +20 dB is a very big boost. Likewise, -3 dB is a small-ish cut, but -60 dB is an extremely deep cut.

All samples are 64-bit floating point. No multiplication necessary in general, but if you're doing FFT then there is some multiplication necessary, yes.

I haven't personally done much with FFT.... that said, I did find this on my hard drive. After a bit of cleanup, what we're left with is an FFT RMS meter. I'll be the first to say that in many ways it's very primitive, and really not exactly that usable.... but maybe it'll give you an idea or two.

desc: FFT slider frequency meter

slider1:-60<-60,30>at 125 Hz
slider2:-60<-60,30>at 250 Hz
slider3:-60<-60,30>at 500 Hz
slider4:-60<-60,30>at 1 kHz
slider5:-60<-60,30>at 2 kHz
slider6:-60<-60,30>at 4 kHz
slider7:-60<-60,30>at 8 kHz
slider8:-60<-60,30>at 16 kHz

// we are setting the FFT size at 512
// this keeps the frequency bins a little wider so it's a fraction bit
// more sensitive, but this does come at the expensive of low-freq
// resolution.

@init 

fftsize=512;
t=0;

freq1 = fftsize*125/srate*2;
freq2 = fftsize*250/srate*2;
freq3 = fftsize*500/srate*2;
freq4 = fftsize*1000/srate*2;
freq5 = fftsize*2000/srate*2;
freq6 = fftsize*4000/srate*2;
freq7 = fftsize*8000/srate*2;
freq8 = fftsize*16000/srate*2;

 bpos=0;
 curblock=7000;
 lastblock=40000;
 window=210000;
 hist=2800000;
 invfsize=1/fftsize;
 hfftsize=fftsize*0.5;
 tmp=0;
 tsc=3.14159/hfftsize;

// this is a window function to help improve FFT spectral performance

 loop(hfftsize,
  window[tmp]=0.42-0.50*cos(tmp*tsc)+0.8*cos(2*tmp*tsc);
  tmp+=1;
  );

 pdc_top_ch=2;
 pdc_bot_ch=0;
 pdc_delay=fftsize;
 ratio2db = 6/log(2);


// I'm hard-coding it to 10 ms
// note that I have to take into account the fact that I'm only
// processing every fftsize samples

 coeff = exp(-1/(10 * 0.001 * srate/fftsize));
 icoeff = 1 - coeff;

function rms(in)
 instance(rms_s)
(
  rms_s = (rms_s * coeff) + (icoeff * in * in);
  sqrt(rms_s);
);
  
@slider

@sample

bpos >= fftsize ? (

 t=curblock;
 curblock=lastblock;
 lastblock=t;

 fft(curblock,fftsize);
 fft_permute(curblock,fftsize);

 slider1 = max(log(s0.rms(curblock[freq1]))*ratio2db,-60);
 slider2 = max(log(s1.rms(curblock[freq2]))*ratio2db,-60);
 slider3 = max(log(s2.rms(curblock[freq3]))*ratio2db,-60);
 slider4 = max(log(s3.rms(curblock[freq4]))*ratio2db,-60);
 slider5 = max(log(s4.rms(curblock[freq5]))*ratio2db,-60);
 slider6 = max(log(s5.rms(curblock[freq6]))*ratio2db,-60);
 slider7 = max(log(s6.rms(curblock[freq7]))*ratio2db,-60);
 slider8 = max(log(s7.rms(curblock[freq8]))*ratio2db,-60);

 i=0;
 loop(hfftsize, 
   i2=fftsize*2-i-2;
   curblock[i] *= invfsize; 
   curblock[i+1] *= invfsize; 
   curblock[i2] *= invfsize; 
   curblock[i2+1] *= invfsize; 
   i+=2; 
   );

 fft_ipermute(curblock,fftsize);
 ifft(curblock,fftsize);
 bpos=0;
 );

 // make sample

 w=window[bpos*0.5];
 iw=1-w;

 os0=spl0;
 os1=spl1;

 spl0=(curblock[bpos]*w + lastblock[fftsize+bpos]*iw);
 spl1=(curblock[bpos+1]*w + lastblock[fftsize+bpos+1]*iw);

 lastblock[bpos]=hist[bpos];
 lastblock[bpos+1]=hist[bpos+1];
 lastblock[fftsize+bpos]=os0;
 lastblock[fftsize+bpos+1]=os1;

 hist[bpos]=os0;
 hist[bpos+1]=os1;
 bpos+=2;

testing this code, btw, is a lot easier if you have a very simple sine oscillator to play with. ahem.

desc:very simple sine oscillator

slider1:500<20,20000,1>Frequency (Hz)
slider2:-15<-60,0,0.1>Volume (dB)


@init

// this is a single pole parameter smoother, so that when you change the frequency slider
// it slides between frequencies instead of making zipper noises

function singlepole_set(ms)
 instance(coeff)
(
  coeff = exp(-1/(ms/1000*srate));
);

function singlepole(in,target)
  instance(coeff)
(
  in*coeff + target*(1-coeff);
);

freq_smooth.singlepole_set(500);
twopi = 2 * $pi;
freq = slider1;

@slider

gain = 2^(slider2/6);

@sample

freq = freq_smooth.singlepole(freq,slider1);
dt = freq/srate;

out = sin(twopi * t) * gain;

spl0 += out;
spl1 += out;

t += dt;

while (t > 1) ( t -= 1; );

[saddle]

Originally Posted by SaulT

ratio = 2^(decibel/6)

This is an approximation of the more proper

ratio = 10^(decibel/20)

It's a pretty decent approximation (within 1% accuracy) until you get to around -50 dB, and even then it's still pretty close (0.0031004 vs 0.0031623). In certain architectures, the first is a little faster. In JS, being interpreted, it doesn't really matter. Probably the only really important thing is to be consistent - either use one or the other, not both.

Anyway, if that needs a little further explanation, decibel is a ratio. Adding 1 dB is the same as multiplying by 1.122, subtracting 6 dB is about the same as dividing by 2. You can think of it as a shorthand, a way to think of audio in an easy linear fashion vs a more conceptually difficult logarithmic one. +1 dB is a small boost. +20 dB is a very big boost. Likewise, -3 dB is a small-ish cut, but -60 dB is an extremely deep cut.

All samples are 64-bit floating point. No multiplication necessary in general, but if you're doing FFT then there is some multiplication necessary, yes.

I haven't personally done much with FFT.... that said, I did find this on my hard drive. After a bit of cleanup, what we're left with is an FFT RMS meter. I'll be the first to say that in many ways it's very primitive, and really not exactly that usable.... but maybe it'll give you an idea or two.

Thanks ashcat_It, and SaulT!

That was throwing me off. I didn't recognize that as a DB conversion. I understand db's, but I only recognize certain math statements associated with it. That helps me a lot.

But as far as working with the samples, if you wanted to apply a -6db volume change to each sample then you would just calculate the ratio as above, and do something like this:

decibel = -6
ratio = 10^(decibel/20)

spl0 = spl0 * ratio
spl1 = spl1 * ratio

and the output would be down -6db.? ( of course with the proper syntax and structure for JS)

Thanks for the FFT code. I'll look through it and try to understand it more.

Saddle

[ashcat_lt]

Originally Posted by saddle

...and the output would be down -6db.? ( of course with the proper syntax and structure for JS?

Yep. I think that part might be a bit obscured in the Volume plug because of the parameter smoothing and whatnot, but it's pretty much what it does. I'm just a hack, so can't help with the FFT part, but I know this much.

[SaulT]

Well, I'll be damned. I started giving a quick summary of the code and before I got anywhere, I realized that the window function was wrong.

window[tmp]=0.42-0.50*cos(tmp*tsc)+0.8*cos(2*tmp*tsc);

should be

window[tmp]=0.42-0.50*cos(tmp*tsc)+0.08*cos(2*tmp*tsc);

It's (supposed to be) a Blackman window....

http://en.wikipedia.org/wiki/Window_...ackman_windows

Okay, that sidetracked me a little. When I have a chance to sit down and think it through I'll see if I can write out a summary of how to FFT or MDCT, at least as best as I understand it.

[SaulT]

Well, I can't say much more than this is an overlap-add FFT method using a Blackman window. It looks like an overlap of 50%. It looks clean when I check it on the spectrograph, so it seems to work correctly.... I don't understand it well enough to design anything better...

[kuus0]

Here's some test code for playing with FFT:

PHP Code:

slider1: 1<0.01,1,.01> t_len (seconds)
slider2: 1<.1,1000,.1> freq (Hz)
slider3: 3.14159265<0,6.283185,0.01> phase (radians)
slider4: 0.5<0,1,0.01> amp

@init
memoff = 1000;

function ar_alloc(nslots) local(ret) (
  ret = memoff;
  memoff += nslots;
  ret;
);

buf_len = 128;
fft_size = buf_len;
in_buf = ar_alloc(buf_len);
tmp_buf = ar_alloc(buf_len);
out_buf = ar_alloc(buf_len);
fft_buf = ar_alloc(buf_len*2);

function add_sine(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += sin(2*$pi*freq*t_len/buf_len*idx + phase) * amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += (((2*idx/(buf_len/(t_len/freq))) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += ((((idx/(buf_len / (2 * t_len*freq))) + (phase/$pi)) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_saw(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  ff = (buf_len / (t_len*freq));
  ph = ff * (.5 * phase/$pi);
  while(idx < buf_len) (
    ar[idx] += ((((idx + ph) % ff) / (ff-1)) * 2 - 1)*amp;
    idx += 1;
  );
);


function analyze_fft(buf, fft_size, t_len) (
  fund_freq = 1/t_len;
  
  hi_bin = 0;
  hi_amp = 0;
  
  idx = 0;
  while(idx < fft_size) (
    (buf[idx*2] > hi_amp) ? (
      hi_bin = idx;
      hi_amp = buf[idx*2];
    );
    idx += 1;
  );
  hi_amp /= fft_size;
  hi_freq = fund_freq * hi_bin;
);

function apply_blackman_window(buf, len, a0,a1,a2,a3) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= a0 
               - a1 * cos(2*$pi*idx/(len-1))
               + a2 * cos(4*$pi*idx/(len-1))
               - a3 * cos(6*$pi*idx/(len-1));
    idx += 1;
  );
);

function apply_blackman_harris(buf, len) (
  apply_blackman_window(buf, len, .35875, .48829, .14128, .01168);
);
  

function do_update() (
  memset(in_buf, 0, buf_len);
  add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq/2, 0, 0.5);
  // add_sine(in_buf, buf_len, t_len, freq, phase, amp);

  memcpy(tmp_buf, in_buf, buf_len);
  apply_blackman_harris(tmp_buf, buf_len);
  
  memset(fft_buf, 0, buf_len*2);
  idx = 0;
  while(idx < buf_len) (
    fft_buf[idx*2] = tmp_buf[idx];
    idx += 1;
  );

  fft(fft_buf, fft_size);
  fft_permute(fft_buf, fft_size);

  analyze_fft(fft_buf, fft_size, t_len); 
  
  memcpy(out_buf, tmp_buf, buf_len); 
);

do_update();

@slider 
t_len = slider1;
freq = slider2;
phase = slider3;
amp = slider4;
do_update();

@gfx 1 1
function setcol(r,g,b) (
  gfx_r = r;
  gfx_g = g;
  gfx_b = b;
);


function draw_wave(ar, len, xo,yo,width,height) (
  setcol(.5,.5,.5);
  gfx_line(xo,yo, xo+width, yo);
  gfx_line(xo,yo + height, xo+width, yo + height);

  setcol(1,1,1);
  pps = width / len;
  spl = 0;
  gfx_x = 0;
  gfx_y = 0;
    
  while(spl < len) (
    // 
    xx = spl * pps + xo;
    yy = height * (ar[spl] / 2 + .5) + yo;
    (gfx_x || gfx_y) ? (
      gfx_lineto(xx,yy);
    ):(
      gfx_x = xx;
      gfx_y = yy;
    );
    (pps > 10) ? (
      gfx_circle(xx,yy,2);
      
    );
    spl += 1;
  );
) ;

setcol(1,1,1);
gfx_x = gfx_y = 0;
gfx_printf("len: %d t_len: %.1fs freq: %.1f phase: %.2f", buf_len, t_len, freq, phase);
gfx_x = 0;
gfx_y += gfx_texth +2;
gfx_printf("fft_size: %d fund_freq: %.2f hi_freq: %.2fHz/%.2f", fft_size, fund_freq, hi_freq, hi_amp);
top = 20;
height = (gfx_h - top) / 2;


setcol(0,1,0);
yhi = top+height;
ylo = top+height*2;
he = height;
wtf= 2;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx, ylo, idx, ylo-(he*(fft_buf[idx*2]/fft_size)*wtf)); 
  idx += 1;
);

draw_wave(in_buf, buf_len, 0,top, gfx_w, height);
draw_wave(out_buf, buf_len, 0,top+height, gfx_w, height); 


I've not used FFT and don't understand it so well, so there may be some fundamental problems with the code. Let me know if I'm doing something wrong.

The upper graph shows the original signal, and lower shows the windowed signal. There's a little FFT spectrogram in the lower left corner.

[saddle]

Originally Posted by SaulT

Well, I'll be damned. I started giving a quick summary of the code and before I got anywhere, I realized that the window function was wrong.

window[tmp]=0.42-0.50*cos(tmp*tsc)+0.8*cos(2*tmp*tsc);

should be

window[tmp]=0.42-0.50*cos(tmp*tsc)+0.08*cos(2*tmp*tsc);

It's (supposed to be) a Blackman window....

http://en.wikipedia.org/wiki/Window_...ackman_windows

Okay, that sidetracked me a little. When I have a chance to sit down and think it through I'll see if I can write out a summary of how to FFT or MDCT, at least as best as I understand it.

Thanks for the personal help. I'll start studying this and the 'window' function as it relates to FFT!

[saddle]

Originally Posted by kuus0

Here's some test code for playing with FFT:

PHP Code:

slider1: 1<0.01,1,.01> t_len (seconds)
slider2: 1<.1,1000,.1> freq (Hz)
slider3: 3.14159265<0,6.283185,0.01> phase (radians)
slider4: 0.5<0,1,0.01> amp

@init
memoff = 1000;

function ar_alloc(nslots) local(ret) (
  ret = memoff;
  memoff += nslots;
  ret;
);

buf_len = 128;
fft_size = buf_len;
in_buf = ar_alloc(buf_len);
tmp_buf = ar_alloc(buf_len);
out_buf = ar_alloc(buf_len);
fft_buf = ar_alloc(buf_len*2);

function add_sine(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += sin(2*$pi*freq*t_len/buf_len*idx + phase) * amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += (((2*idx/(buf_len/(t_len/freq))) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += ((((idx/(buf_len / (2 * t_len*freq))) + (phase/$pi)) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_saw(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  ff = (buf_len / (t_len*freq));
  ph = ff * (.5 * phase/$pi);
  while(idx < buf_len) (
    ar[idx] += ((((idx + ph) % ff) / (ff-1)) * 2 - 1)*amp;
    idx += 1;
  );
);


function analyze_fft(buf, fft_size, t_len) (
  fund_freq = 1/t_len;
  
  hi_bin = 0;
  hi_amp = 0;
  
  idx = 0;
  while(idx < fft_size) (
    (buf[idx*2] > hi_amp) ? (
      hi_bin = idx;
      hi_amp = buf[idx*2];
    );
    idx += 1;
  );
  hi_amp /= fft_size;
  hi_freq = fund_freq * hi_bin;
);

function apply_blackman_window(buf, len, a0,a1,a2,a3) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= a0 
               - a1 * cos(2*$pi*idx/(len-1))
               + a2 * cos(4*$pi*idx/(len-1))
               - a3 * cos(6*$pi*idx/(len-1));
    idx += 1;
  );
);

function apply_blackman_harris(buf, len) (
  apply_blackman_window(buf, len, .35875, .48829, .14128, .01168);
);
  

function do_update() (
  memset(in_buf, 0, buf_len);
  add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq/2, 0, 0.5);
  // add_sine(in_buf, buf_len, t_len, freq, phase, amp);

  memcpy(tmp_buf, in_buf, buf_len);
  apply_blackman_harris(tmp_buf, buf_len);
  
  memset(fft_buf, 0, buf_len*2);
  idx = 0;
  while(idx < buf_len) (
    fft_buf[idx*2] = tmp_buf[idx];
    idx += 1;
  );

  fft(fft_buf, fft_size);
  fft_permute(fft_buf, fft_size);

  analyze_fft(fft_buf, fft_size, t_len); 
  
  memcpy(out_buf, tmp_buf, buf_len); 
);

do_update();

@slider 
t_len = slider1;
freq = slider2;
phase = slider3;
amp = slider4;
do_update();

@gfx 1 1
function setcol(r,g,b) (
  gfx_r = r;
  gfx_g = g;
  gfx_b = b;
);


function draw_wave(ar, len, xo,yo,width,height) (
  setcol(.5,.5,.5);
  gfx_line(xo,yo, xo+width, yo);
  gfx_line(xo,yo + height, xo+width, yo + height);

  setcol(1,1,1);
  pps = width / len;
  spl = 0;
  gfx_x = 0;
  gfx_y = 0;
    
  while(spl < len) (
    // 
    xx = spl * pps + xo;
    yy = height * (ar[spl] / 2 + .5) + yo;
    (gfx_x || gfx_y) ? (
      gfx_lineto(xx,yy);
    ):(
      gfx_x = xx;
      gfx_y = yy;
    );
    (pps > 10) ? (
      gfx_circle(xx,yy,2);
      
    );
    spl += 1;
  );
) ;

setcol(1,1,1);
gfx_x = gfx_y = 0;
gfx_printf("len: %d t_len: %.1fs freq: %.1f phase: %.2f", buf_len, t_len, freq, phase);
gfx_x = 0;
gfx_y += gfx_texth +2;
gfx_printf("fft_size: %d fund_freq: %.2f hi_freq: %.2fHz/%.2f", fft_size, fund_freq, hi_freq, hi_amp);
top = 20;
height = (gfx_h - top) / 2;


setcol(0,1,0);
yhi = top+height;
ylo = top+height*2;
he = height;
wtf= 2;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx, ylo, idx, ylo-(he*(fft_buf[idx*2]/fft_size)*wtf)); 
  idx += 1;
);

draw_wave(in_buf, buf_len, 0,top, gfx_w, height);
draw_wave(out_buf, buf_len, 0,top+height, gfx_w, height); 


I've not used FFT and don't understand it so well, so there may be some fundamental problems with the code. Let me know if I'm doing something wrong.

The upper graph shows the original signal, and lower shows the windowed signal. There's a little FFT spectrogram in the lower left corner.

WOW! Thank you! I'll load this up and start trying it out.

Where I'm lost so to speak right now, is that I understand how the FFT transfers from a Time to a Frequency based view. And that the information is contained in 'frequency bins' that vary in size based on the sampleblock size and sample rate.

It's the process now of trying to understand what I can do with the bare PCM samples in terms of levels, filters, etc. and the math to do that properly, and then the next step, understanding when and how to use the data returned in the 'frequency bins.'

A lot to learn. Thank you all for your help!
Saddle

[kuus0]

There's a bug in my code, where the signals are drawn flipped.

Add a minus sign as shown below in red:

yy = height * (-ar[spl] / 2 + .5) + yo;


Here's a second version. I'm trying to generate the original signal from the FFT data, but I'm not getting very promising results. There's something very fishy going on. If you play with the input signals phase (3rd slider) the FFT result changes like crazy.

The recovered signal is in red.

PHP Code:

slider1: 1<0.01,1,.01> t_len (seconds)
slider2: 10<.1,1000,.1> freq (Hz)
slider3: 1.570796<0,6.283185,0.01> phase (radians)
slider4: 0.9<0,1,0.01> amp
slider5: 2<0, 5,1{32,64,128,256,512,1024}> fft_size
//3.14159265
@init
memoff = 1000;

function ar_alloc(nslots) local(ret) (
  ret = memoff;
  memoff += nslots;
  ret;
);

max_buf_len = 2^10;
in_buf = ar_alloc (max_buf_len);
tmp_buf = ar_alloc(max_buf_len);
out_buf = ar_alloc(max_buf_len);
fft_buf = ar_alloc(max_buf_len*2);

function add_sine(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += sin(2*$pi*freq*t_len/buf_len*idx + phase) * amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += (((2*idx/(buf_len/(t_len/freq))) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += ((((idx/(buf_len / (2 * t_len*freq))) + (phase/$pi)) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_saw(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  ff = (buf_len / (t_len*freq));
  ph = ff * (.5 * phase/$pi);
  while(idx < buf_len) (
    ar[idx] += ((((idx + ph) % ff) / (ff-1)) * 2 - 1)*amp;
    idx += 1;
  );
);


function analyze_fft(buf, fft_size, t_len) (
  fund_freq = 1/t_len;
  
  hi_bin = 0;
  hi_amp = 0;
  
  idx = 0;
  while(idx < fft_size) (
    (buf[idx*2] > hi_amp) ? (
      hi_bin = idx;
      hi_amp = buf[idx*2];
    );
    idx += 1;
  );
  hi_amp /= fft_size;
  hi_freq = fund_freq * hi_bin;
);

function apply_blackman_window(buf, len, a0,a1,a2,a3) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= a0 
               - a1 * cos(2*$pi*idx/(len-1))
               + a2 * cos(4*$pi*idx/(len-1))
               - a3 * cos(6*$pi*idx/(len-1));
    idx += 1;
  );
);

function compute_from_fft(fft_buf, fft_size, t_len, x_s) local(idx)  (
  v = 0;
  fund_f = 1/t_len;
  idx = 0; 
  while(idx <=  floor(fft_size / 2)) (
    v += sin(x_s * fund_f * (idx+1) * (2*$pi) + (fft_buf[idx*2+1]/fft_size)) * (fft_buf[idx*2] / fft_size);
    idx += 1; 
  );
  v;
);

function apply_blackman_harris(buf, len) (
  apply_blackman_window(buf, len, .35875, .48829, .14128, .01168);
);
  

function do_update() (
  memset(in_buf, 0, buf_len);
  //add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq/2, 0, 0.5);
  add_sine(in_buf, buf_len, t_len, freq, phase, amp);

  memcpy(tmp_buf, in_buf, buf_len);
  apply_blackman_harris(tmp_buf, buf_len);
  
  memset(fft_buf, 0, buf_len*2);
  idx = 0;
  while(idx < buf_len) (
    fft_buf[idx*2] = tmp_buf[idx];
    idx += 1;
  );

  fft(fft_buf, fft_size);
  fft_permute(fft_buf, fft_size);

  analyze_fft(fft_buf, fft_size, t_len); 
  
  memcpy(out_buf, tmp_buf, buf_len); 
);

do_update();

@slider 
t_len = slider1;
freq = slider2;
phase = slider3;
amp = slider4;

buf_len = 2^(slider5+5);
fft_size = 2^(slider5+5);
do_update();

@gfx 1 1
function setcol(r,g,b) (
  gfx_r = r;
  gfx_g = g;
  gfx_b = b;
);


function draw_wave(ar, len, xo,yo,width,height) (
  setcol(.5,.5,.5);
  gfx_line(xo,yo, xo+width, yo);
  gfx_line(xo,yo + height, xo+width, yo + height);

  setcol(1,1,1);
  pps = width / len;
  spl = 0;
  gfx_x = 0;
  gfx_y = 0;
    
  while(spl < len) (
    // 
    xx = spl * pps + xo;
    yy = height * (-ar[spl] / 2 + .5) + yo;
    (gfx_x || gfx_y) ? (
      gfx_lineto(xx,yy);
    ):(
      gfx_x = xx;
      gfx_y = yy;
    );
    (pps > 10) ? (
      gfx_circle(xx,yy,2);
      
    );
    spl += 1;
  );
);

function draw_from_fft(fft_buf, fft_size, xo,yo,width,height) (
  gfx_x = xo;
  gfx_y = (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width*t_len))) / 2 * height + yo;
  while(gfx_x < (width+xo)) (
    //
    gfx_lineto(gfx_x, (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width*t_len))) / 2 * height + yo);
    gfx_x += 1; 
  );
);

setcol(1,1,1);
gfx_x = gfx_y = 0;
gfx_printf("len: %d t_len: %.1fs freq: %.1f phase: %.2f", buf_len, t_len, freq, phase);
gfx_x = 0;
gfx_y += gfx_texth +2;
gfx_printf("fft_size: %d fund_freq: %.2f hi_freq: %.2fHz/%.2f", fft_size, fund_freq, hi_freq, hi_amp);
top = 20;
height = (gfx_h - top) / 2;


setcol(0,1,0);
yhi = top+height;
ylo = top+height*2;
he = height;
wtf= 2;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx, ylo, idx, ylo-(he*(fft_buf[idx*2]/fft_size)*wtf)); 
  idx += 1;
);

draw_wave(in_buf, buf_len, 0,top, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top, gfx_w-6, height);

draw_wave(out_buf, buf_len, 0,top+height, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top+height, gfx_w-6, height); 


[witti]

I was reading this thread with interest. Nice scripts, folks !

Here is a challenge for you: A center canceller plugin !

Unfortunately i'm not able to embed the FFT stuff. Here is the raw code.
It's from the musicdsp site: http://www.musicdsp.org/showArchiveC...?ArchiveID=163

slider1:1<0,1,1>Left
slider2:0<0,1,1>Center
slider3:1<0,1,1>Right

@init

@slider
AmpL=slider1;
AmpC=slider2;
AmpR=slider3;

@sample
ModL = spl0 [i] * spl0 [i] + spl0 [j] * spl0 [j];
ModR = spl1 [i] * spl1 [i] + spl1 [j] * spl1 [j];

// min on complex numbers
( ModL > ModR ) ?
(
RealC = spl1 [i];
ImagC = spl1 [j];
):(
RealC = spl0 [i];
ImagC = spl0 [j];
);

// phase correction...
u = abs(atan2(spl0 [j], spl0 [i]) - atan2(spl1 [j], spl1 [i])) / $pi;

( u >= 1 ) ? u = 2 - u;

u = pow(1 - u*u*u, 24);

RealC *= u;
ImagC *= u;

// center extraction...
RealL = spl0 [i] - RealC;
ImagL = spl0 [j] - ImagC;

RealR = spl1 [i] - RealC;
ImagR = spl1 [j] - ImagC;

// You can do some treatments here...

spl0 [i] = RealL * AmpL + RealC * AmpC;
spl0 [j] = ImagL * AmpL + ImagC * AmpC;

spl1 [i] = RealR * AmpR + RealC * AmpC;
spl1 [j] = ImagR * AmpR + ImagC * AmpC;

Wouldn't it be nice to have a center canceller for reaper ?

[saddle]

Originally Posted by kuus0

There's a bug in my code, where the signals are drawn flipped.

Add a minus sign as shown below in red:

yy = height * (-ar[spl] / 2 + .5) + yo;


Here's a second version. I'm trying to generate the original signal from the FFT data, but I'm not getting very promising results. There's something very fishy going on. If you play with the input signals phase (3rd slider) the FFT result changes like crazy.

The recovered signal is in red.

PHP Code:

slider1: 1<0.01,1,.01> t_len (seconds)
slider2: 10<.1,1000,.1> freq (Hz)
slider3: 1.570796<0,6.283185,0.01> phase (radians)
slider4: 0.9<0,1,0.01> amp
slider5: 2<0, 5,1{32,64,128,256,512,1024}> fft_size
//3.14159265
@init
memoff = 1000;

function ar_alloc(nslots) local(ret) (
  ret = memoff;
  memoff += nslots;
  ret;
);

max_buf_len = 2^10;
in_buf = ar_alloc (max_buf_len);
tmp_buf = ar_alloc(max_buf_len);
out_buf = ar_alloc(max_buf_len);
fft_buf = ar_alloc(max_buf_len*2);

function add_sine(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += sin(2*$pi*freq*t_len/buf_len*idx + phase) * amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += (((2*idx/(buf_len/(t_len/freq))) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += ((((idx/(buf_len / (2 * t_len*freq))) + (phase/$pi)) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_saw(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  ff = (buf_len / (t_len*freq));
  ph = ff * (.5 * phase/$pi);
  while(idx < buf_len) (
    ar[idx] += ((((idx + ph) % ff) / (ff-1)) * 2 - 1)*amp;
    idx += 1;
  );
);


function analyze_fft(buf, fft_size, t_len) (
  fund_freq = 1/t_len;
  
  hi_bin = 0;
  hi_amp = 0;
  
  idx = 0;
  while(idx < fft_size) (
    (buf[idx*2] > hi_amp) ? (
      hi_bin = idx;
      hi_amp = buf[idx*2];
    );
    idx += 1;
  );
  hi_amp /= fft_size;
  hi_freq = fund_freq * hi_bin;
);

function apply_blackman_window(buf, len, a0,a1,a2,a3) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= a0 
               - a1 * cos(2*$pi*idx/(len-1))
               + a2 * cos(4*$pi*idx/(len-1))
               - a3 * cos(6*$pi*idx/(len-1));
    idx += 1;
  );
);

function compute_from_fft(fft_buf, fft_size, t_len, x_s) local(idx)  (
  v = 0;
  fund_f = 1/t_len;
  idx = 0; 
  while(idx <=  floor(fft_size / 2)) (
    v += sin(x_s * fund_f * (idx+1) * (2*$pi) + (fft_buf[idx*2+1]/fft_size)) * (fft_buf[idx*2] / fft_size);
    idx += 1; 
  );
  v;
);

function apply_blackman_harris(buf, len) (
  apply_blackman_window(buf, len, .35875, .48829, .14128, .01168);
);
  

function do_update() (
  memset(in_buf, 0, buf_len);
  //add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq/2, 0, 0.5);
  add_sine(in_buf, buf_len, t_len, freq, phase, amp);

  memcpy(tmp_buf, in_buf, buf_len);
  apply_blackman_harris(tmp_buf, buf_len);
  
  memset(fft_buf, 0, buf_len*2);
  idx = 0;
  while(idx < buf_len) (
    fft_buf[idx*2] = tmp_buf[idx];
    idx += 1;
  );

  fft(fft_buf, fft_size);
  fft_permute(fft_buf, fft_size);

  analyze_fft(fft_buf, fft_size, t_len); 
  
  memcpy(out_buf, tmp_buf, buf_len); 
);

do_update();

@slider 
t_len = slider1;
freq = slider2;
phase = slider3;
amp = slider4;

buf_len = 2^(slider5+5);
fft_size = 2^(slider5+5);
do_update();

@gfx 1 1
function setcol(r,g,b) (
  gfx_r = r;
  gfx_g = g;
  gfx_b = b;
);


function draw_wave(ar, len, xo,yo,width,height) (
  setcol(.5,.5,.5);
  gfx_line(xo,yo, xo+width, yo);
  gfx_line(xo,yo + height, xo+width, yo + height);

  setcol(1,1,1);
  pps = width / len;
  spl = 0;
  gfx_x = 0;
  gfx_y = 0;
    
  while(spl < len) (
    // 
    xx = spl * pps + xo;
    yy = height * (-ar[spl] / 2 + .5) + yo;
    (gfx_x || gfx_y) ? (
      gfx_lineto(xx,yy);
    ):(
      gfx_x = xx;
      gfx_y = yy;
    );
    (pps > 10) ? (
      gfx_circle(xx,yy,2);
      
    );
    spl += 1;
  );
);

function draw_from_fft(fft_buf, fft_size, xo,yo,width,height) (
  gfx_x = xo;
  gfx_y = (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width*t_len))) / 2 * height + yo;
  while(gfx_x < (width+xo)) (
    //
    gfx_lineto(gfx_x, (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width*t_len))) / 2 * height + yo);
    gfx_x += 1; 
  );
);

setcol(1,1,1);
gfx_x = gfx_y = 0;
gfx_printf("len: %d t_len: %.1fs freq: %.1f phase: %.2f", buf_len, t_len, freq, phase);
gfx_x = 0;
gfx_y += gfx_texth +2;
gfx_printf("fft_size: %d fund_freq: %.2f hi_freq: %.2fHz/%.2f", fft_size, fund_freq, hi_freq, hi_amp);
top = 20;
height = (gfx_h - top) / 2;


setcol(0,1,0);
yhi = top+height;
ylo = top+height*2;
he = height;
wtf= 2;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx, ylo, idx, ylo-(he*(fft_buf[idx*2]/fft_size)*wtf)); 
  idx += 1;
);

draw_wave(in_buf, buf_len, 0,top, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top, gfx_w-6, height);

draw_wave(out_buf, buf_len, 0,top+height, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top+height, gfx_w-6, height); 


Hey! That's some cool code! It works, and now I have to disect it and figure out what it's doing. Something I've really been wondering. How/when to use the data from the FFT transform in the frequency domain.

Thank you!

[kuus0]

I played with this a bit more and now the reverse from the FFT seems to work, yay!

I used Eq. 2 from here:
https://en.wikipedia.org/wiki/Discre...rier_transform

It's a sum of products of complex numbers. Then I used Euler's formula on the exponent thing. Then, because we know the original data is real only, the we can simplify the complex number product (throw out the imaginary part, it should be 0). This is all happens in compute_from_fft.

PHP Code:

slider1: 1<0.01,1,.01> t_len (seconds)
slider2: 10<.1,1000,.1> freq (Hz)
slider3: 1.570796<0,6.283185,0.01> phase (radians)
slider4: 0.9<0,1,0.01> amp
slider5: 2<0, 5,1{32,64,128,256,512,1024}> fft_size
slider6: 0<0, 3,1{sine,square,saw}> osc shape
slider7: 0<0, 2,1{blackman harris, hann, rectangular}> window shape
//3.14159265
@init
memoff = 1000;

function ar_alloc(nslots) local(ret) (
  ret = memoff;
  memoff += nslots;
  ret;
);

max_buf_len = 2^10;
in_buf = ar_alloc (max_buf_len);
tmp_buf = ar_alloc(max_buf_len);
out_buf = ar_alloc(max_buf_len);
fft_buf = ar_alloc(max_buf_len*2);

function add_sine(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += sin(2*$pi*freq*t_len/buf_len*idx + phase) * amp;
    idx += 1;
  );
);

function add_square(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  while(idx < buf_len) (
    ar[idx] += ((((idx/(buf_len / (2 * t_len*freq))) + (phase/$pi)) % 2) * 2 - 1)*amp;
    idx += 1;
  );
);

function add_saw(ar, buf_len, t_len, freq, phase, amp) (
  idx = 0;
  ff = (buf_len / (t_len*freq));
  ph = ff * (.5 * phase/$pi);
  while(idx < buf_len) (
    ar[idx] += ((((idx + ph) % ff) / (ff-1)) * 2 - 1)*amp;
    idx += 1;
  );
);


function analyze_fft(buf, fft_size, t_len) (
  fund_freq = 1/t_len;
  
  hi_bin = 0;
  hi_amp = 0;
  
  idx = 0;
  while(idx < fft_size) (
    (buf[idx*2] > hi_amp) ? (
      hi_bin = idx;
      hi_amp = buf[idx*2];
    );
    idx += 1;
  );
  hi_amp /= fft_size;
  hi_freq = fund_freq * hi_bin;
);

function apply_blackman_window(buf, len, a0,a1,a2,a3) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= a0 
               - a1 * cos(2*$pi*idx/(len-1))
               + a2 * cos(4*$pi*idx/(len-1))
               - a3 * cos(6*$pi*idx/(len-1));
    idx += 1;
  );
);

function apply_hann_window(buf, len) (
  idx = 0;
  while(idx < len) (
    //
    buf[idx] *= .5 * (1 - cos((2*$pi*idx/(len-1)))); 
    idx += 1;
  );
);


function compute_from_fft(fft_buf, fft_size, t_len, x_rel) local(idx)  (
  v = 0;
  idx = 0;
  scalefactor = 2 * (1/fft_size);
  while(idx < floor(fft_size / 2)) (
    tt = x_rel * idx * (2*$pi);
    v += cos(tt) * fft_buf[idx*2] - fft_buf[idx*2+1] * sin(tt);
    idx += 1; 
  );
  v * scalefactor;
);

function apply_blackman_harris(buf, len) (
  apply_blackman_window(buf, len, .35875, .48829, .14128, .01168);
);
  

function do_update() (
  memset(in_buf, 0, buf_len);
  (slider6 === 0) ? add_sine(in_buf, buf_len, t_len, freq, phase, amp);
  (slider6 === 1) ? add_square(in_buf, buf_len, t_len, freq, phase, amp);
  (slider6 === 2) ? add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  //add_saw(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq, phase, amp);
  // add_square(in_buf, buf_len, t_len, freq/2, 0, 0.5);
  //

  memcpy(tmp_buf, in_buf, buf_len);
  (slider7 === 0) ? apply_blackman_harris(tmp_buf, buf_len);
  (slider7 === 1) ? apply_hann_window(tmp_buf, buf_len);
  (slider7 === 2) ? 1; // rectangular 
  
  memset(fft_buf, 0, buf_len*2);
  idx = 0;
  while(idx < buf_len) (
    fft_buf[idx*2] = tmp_buf[idx];
    idx += 1;
  );

  fft(fft_buf, fft_size);
  fft_permute(fft_buf, fft_size);

  analyze_fft(fft_buf, fft_size, t_len); 
  
  memcpy(out_buf, tmp_buf, buf_len); 
);

do_update();

@slider 
t_len = slider1;
freq = slider2;
phase = slider3;
amp = slider4;

buf_len = 2^(slider5+5);
fft_size = buf_len;
do_update();

@gfx 1 1
function setcol(r,g,b) (
  gfx_r = r;
  gfx_g = g;
  gfx_b = b;
);


function draw_wave(ar, len, xo,yo,width,height) (
  setcol(.5,.5,.5);
  gfx_line(xo,yo, xo+width, yo);
  gfx_line(xo,yo + height, xo+width, yo + height);

  setcol(1,1,1);
  pps = width / len;
  spl = 0;
  gfx_x = 0;
  gfx_y = 0;
    
  while(spl < len) (
    // 
    xx = spl * pps + xo;
    yy = height * (-ar[spl] / 2 + .5) + yo;
    (gfx_x || gfx_y) ? (
      gfx_lineto(xx,yy);
    ):(
      gfx_x = xx;
      gfx_y = yy;
    );
    (pps > 10) ? (
      gfx_circle(xx,yy,2);
      
    );
    spl += 1;
  );
);

function draw_from_fft(fft_buf, fft_size, xo,yo,width,height) (
  gfx_x = xo;
  gfx_y = (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width*t_len))) / 2 * height + yo;
  while(gfx_x < (width+xo)) (
    //
    yy = (1+compute_from_fft(fft_buf, fft_size, t_len, (gfx_x/width))) / 2;
    gfx_lineto(gfx_x, -yy * height + yo+height); 
    gfx_x += 1; 
  );
);

setcol(1,1,1);
gfx_x = gfx_y = 0;
gfx_printf("len: %d t_len: %.1fs freq: %.1f phase: %.2f", buf_len, t_len, freq, phase);
gfx_x = 0;
gfx_y += gfx_texth +2;
nyquist = 1/(t_len/fft_size)/2;
gfx_printf("fft_size: %d fund_freq: %.2fHz nyquist: %.2fHz hi_freq: %.2fHz/%.2f", 
           fft_size,
           fund_freq,
           nyquist,
           hi_freq,
           hi_amp);
top = 20;
height = (gfx_h - top) / 2;


setcol(0,1,0);
yhi = top+height;
ylo = top+height*2;
he = height;
wtf= 2;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx, ylo, idx, ylo-(he*(fft_buf[idx*2]/fft_size)*wtf)); 
  idx += 1;
);

setcol(0,1,1);
wtf= 1;
xo = gfx_x;
idx = 0;
while(idx < fft_size) (
  gfx_line(idx + xo, ylo, idx + xo, ylo-(he*(fft_buf[idx*2+1]/fft_size)*wtf)); 
  idx += 1;
);

draw_wave(in_buf, buf_len, 0,top, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top, gfx_w-6, height);

draw_wave(out_buf, buf_len, 0,top+height, gfx_w, height);
setcol(1,0,0);
draw_from_fft(fft_buf, fft_size, 3 ,top+height, gfx_w-6, height); 


[kuus0]

Originally Posted by witti

I was reading this thread with interest. Nice scripts, folks !

Here is a challenge for you: A center canceller plugin !

Unfortunately i'm not able to embed the FFT stuff. Here is the raw code.
It's from the musicdsp site: http://www.musicdsp.org/showArchiveC...?ArchiveID=163
[snip]
Wouldn't it be nice to have a center canceller for reaper ?

Interesting! How does this differ from removing the center using mid/side processing?

Do you have sound samples?

[SaulT]

Okay, so here is my now tested and moderately unoptimized stab at this plugin. Finally got it working, and as far as I can tell I'm implementing it as shown in the source, but I'm not really impressed with the results. Unfortunately I don't know the underlying math that he's using well enough to really tweak it.

Feel free to play with it, though, YMMV.

desc:center cancel plugin
slider1:1<0,1,0.01>cancellation (0 = none, 1 = full)
slider2:3<0,4,1{512,1024,2048,4096,8192}>FFT size

@init

oslider2 = -1;

@slider

fSide = cos(slider1/2*$pi);
fMid = cos((1-slider1)/2*$pi);

oslider2 != slider2 ? (
 fftsize = 2^(9+slider2);
 bpos=0;
 curblock=7000;
 lastblock=40000;
 window=210000;
 hist=2800000;
 invfsize=1/fftsize;
 hfftsize=fftsize*0.5;
 tmp=0;
 tsc=$pi/hfftsize;

 loop(hfftsize,
  window[tmp]=0.42-0.50*cos(tmp*tsc)+0.08*cos(2*tmp*tsc);
  tmp+=1;
  );

 pdc_top_ch=2;
 pdc_bot_ch=0;
 pdc_delay=fftsize;
 oslider2 = slider2;
 );

@sample

bpos >= fftsize ? (

 t=curblock;
 curblock=lastblock;
 lastblock=t;

 fft(curblock,fftsize);
 fft_permute(curblock,fftsize);

 i=0;
 loop(hfftsize, 
   i2=fftsize*2-i-2;
   curblock[i] *= invfsize; 
   curblock[i+1] *= invfsize; 
   curblock[i2] *= invfsize; 
   curblock[i2+1] *= invfsize; 
   i+=2; 
   );

 fmodL = 0; fmodR = 0;
 fRealL = 0; fRealR = 0; fRealC = 0;
 fImagL = 0; fImagR = 0; fImagC = 0;
 u = 0;


 i0 = 0; i1 = 1;
 j0 = hfftsize; j1 = hfftsize+1;
 loop(fftsize,


   fModL = curblock[i0]*curblock[i0] + curblock[j0]*curblock[j0];
   fModR = curblock[i1]*curblock[i1] + curblock[j1]*curblock[j1];

   (fModL > fModR) ? (
     fRealC = curblock[i1];
     fImagC = curblock[j1];
     ) : (
     fRealC = curblock[i0];
     fImagC = curblock[j0];
     );

   u = abs(atan2(curblock[j0],curblock[i0]) - atan2(curblock[j1],curblock[i1]))/$pi;

   debug1 = atan2(curblock[j0],curblock[i0]);
   debug2 = atan2(curblock[j1],curblock[i1]);

   (u >= 1) ? u = 2 - u;
   u = pow(1-u*u*u,24);

   fRealC = fRealC * u;
   fImagC = fImagC * u;

   fRealL = (curblock[i0] - fRealC) * fSide;
   fImagL = (curblock[j0] - fImagC) * fSide;

   fRealR = (curblock[i1] - fRealC) * fSide;
   fImagR = (curblock[j1] - fImagC) * fSide;

   curblock[i0] = fRealL + fRealC*fMid;
   curblock[j0] = fImagL + fImagC*fMid;

   curblock[i1] = fRealR + fRealC*fMid;
   curblock[j1] = fImagR + fImagC*fMid;

   i0 += 2; j0 += 2;
   i1 += 2; j1 += 2;
   );

 fft_ipermute(curblock,fftsize);
 ifft(curblock,fftsize);
 bpos=0;
 );

 // make sample

 w=window[bpos*0.5];
 iw=1-w;

 os0=spl0;
 os1=spl1;

 spl0=(curblock[bpos]*w + lastblock[fftsize+bpos]*iw);
 spl1=(curblock[bpos+1]*w + lastblock[fftsize+bpos+1]*iw);

 lastblock[bpos]=hist[bpos];
 lastblock[bpos+1]=hist[bpos+1];
 lastblock[fftsize+bpos]=os0;
 lastblock[fftsize+bpos+1]=os1;

 hist[bpos]=os0;
 hist[bpos+1]=os1;
 bpos+=2;

[saddle]

So I think that i understand some...

the fft() function does the transform from time to frequency domain... ?

what does the fft_permute() function do next?

Trying to work through the code.

Thanks!

[saddle]

Originally Posted by saddle

So I think that i understand some...

the fft() function does the transform from time to frequency domain... ?

what does the fft_permute() function do next?

Trying to work through the code.

Thanks!

I think I figured out that the fft_permute flips the frequency order of the bins from high->low to low->high the way we expect it.

Found it in a thread from 2008.

[witti]

Thanks, SaulT. Tested you plugin. (Not working yet.) Had a look at the code. I almost tried the same to implement the fft stuff, borrowing code snipppets from other js fft plugins, but with no luck. At some point i always got stuck.

Sorry Saddle, don't want to hijack your thread. Thought it would be interesting for you, too.

witti

[saddle]

Originally Posted by witti

Thanks, SaulT. Tested you plugin. (Not working yet.) Had a look at the code. I almost tried the same to implement the fft stuff, borrowing code snipppets from other js fft plugins, but with no luck. At some point i always got stuck.

Sorry Saddle, don't want to hijack your thread. Thought it would be interesting for you, too.

witti

No problem! I'm enjoying the mental exercise

[saddle]

I downloaded Sage Math and Virtual Box and installed it this morning. I want to try using it to plot some of the math to visualize what is being done to the signal in the Frequency domain.

It's free. Installs very easy.

www.sagemath.org

Of course, it has a learning curve. Maybe I can plot that.

Found this resource:
http://www.arachnoid.com/sage/fourier.html

with an FFT Sage Math resource file. Just load it into Sage and you have all the math that is included in the article to work with.

[SaulT]

Okay, after quite a bit of monkeying around and fixing the things I got wrong in the above code, I have a functional plugin. As far as I can tell I am implementing the algorithm correctly, but I'm not impressed with the results.

Code is updated above.

[witti]

Thanks SaulT ! I think methods like these will always produce 'smearings' and what not. The result could be better, i agree. At least i've learned something new and your code showed me where i made my mistakes.
Have downloaded the source code of the dsp_centercut (winamp,foobar ?).
I doubt that i can port this to js, but i will have a look and try it...

[kuus0]

Originally Posted by SaulT

Okay, after quite a bit of monkeying around and fixing the things I got wrong in the above code, I have a functional plugin. As far as I can tell I am implementing the algorithm correctly, but I'm not impressed with the results.

Code is updated above.

I think you should be doing fft twice, once for each channel.

[SaulT]

I think you should be doing fft twice, once for each channel.

I would normally have done that from the beginning, but every bundled FFT example uses interleaved stereo, and it works - stereo information is retained from before/after the FFT. I've been trying to understand the following, for instance, but about all I can glean from it is that the theory is solid and backs up what my ears are telling me.

http://www.fftguru.com/fftguru.com.tutorial2.pdf

I downloaded Sage Math and Virtual Box and installed it this morning. I want to try using it to plot some of the math to visualize what is being done to the signal in the Frequency domain.

Plot everything. I use something called Flotr to graph just about every function that I work on because sometimes you need to see the curve to understand something about it. When I'm doing approximations, I will normally pull it out and start graphing like crazy. I've now thrown away one of my tanh() approximations because I started using it before I graphed both the positive and negative sides of the curve, and it's also how I realized that I had the wrong clamp values on my function above.

Which exact tool you use isn't quite as important as long as you have a sandbox environment to evaluate the function.

[kuus0]

Originally Posted by SaulT

I would normally have done that from the beginning, but every bundled FFT example uses interleaved stereo, and it works - stereo information is retained from before/after the FFT. I've been trying to understand the following, for instance, but about all I can glean from it is that the theory is solid and backs up what my ears are telling me.

http://www.fftguru.com/fftguru.com.tutorial2.pdf

Wow, I'm ignorant, as usual I had a quick look and my head started to hurt. Have to try more seriously later. Thanks for the link!

I tried the center canceller implementation. I used mp3 audio and it sounded like it made everything center mono with really heavy (mp3-) compression like artifacts.

There's one center canceller JS maybe it was LOSER's. Should these two sound similar?