How do you do linear phase filters?

[IXix]

I'd like to combine linear phase hi and lo pass filters in a single JS. Ideally I'd want non-resonant filters with params for cutoff and steepness, if such a thing is possible (am I asking for the moon on a stick? my DSP knowledge is pretty thin.)

Can anybody post some code, explain how it's done or point me to some suitable resource? Looking around the net it seems I need to use FIR which is easy enough but how do I work out the necessary coefficients to get the desired frequency response?

Thanks in advance.

[bang]

hi IXix,

you might want to look at bangzero_wsfilter.jsfx-inc, which is included in wild wave fu. it is a port of some windowed sync filter code from musicdsp.org, which was the simplest FIR filter i could find. it does have logic to generate coefficents for bandpass (hi+lowpass) filters. cutoff freq is a simple parameter specifying a fraction of the sample frequency. steepness is less clear: it is determined by the number of coefficents in the filter, which does not map to dB/octave in an obvious way. empirically, each 24 coefficients give about 12 dB/oct, iirc. the biggest issue i had was that computing the filter from input/coefficients in array mem was pretty slow. so you will also find some code that optimizes the filter using instance variables. the code is not well commented, esp the optimized bits. but i could perhaps fix that. and of course i'm happy to answer questions. after all, we dsp hackers need to stick together. those dsp math folks are a tough crowd. :^)

enjoy! /dan

[IXix]

Thanks Dan, I'll check that out. With a good tail-wind I might even be able to understand some of it. Brace yourself for a brace of questions.

[IXix]

The windowed sinc code looks like exactly what I need but I'm not sure how to actually use it. Any chance you could point me in the right direction?

[bang]

Originally Posted by IXix

The windowed sinc code looks like exactly what I need but I'm not sure how to actually use it. Any chance you could point me in the right direction?

sure. give me a day or two and i'll put together a demo jsfx.

[IXix]

Originally Posted by bang

sure. give me a day or two and i'll put together a demo jsfx.

You're a star. Thanks!

[Tale]

Meanwhile you will find some IMHO useful information here:
http://www.labbookpages.co.uk/audio/firWindowing.html

[IXix]

Originally Posted by Tale

Meanwhile you will find some IMHO useful information here:
http://www.labbookpages.co.uk/audio/firWindowing.html

Thanks Tale, that is very helpful. If I read it another twenty or so times I might even understand it.

[saul_sabia]

If it helps, there are examples of windowed-sinc filters at

http://www.dspguide.com/ch16/4.htm

In this case you'll need to know enough Basic to translate to JS. Fc is the cut-off frequency, and is normalized as a percentage of the sampling rate, from 0 to 0.5 . I won't claim that it's the best method, nor that it's optimized, but it should give an example to go from, to start at least.

The hardest part for me to remember as I'm porting code over is that JS only has one array, and having to tweak the code to where the variable is the pointer and not a separate array.

Does that help at all?

[IXix]

Originally Posted by saul_sabia

If it helps, there are examples of windowed-sinc filters at

http://www.dspguide.com/ch16/4.htm
...

Yes, that's helpful too. Thanks.

[Tale]

Originally Posted by IXix

Thanks Tale, that is very helpful. If I read it another twenty or so times I might even understand it.

Well, actually it isn't that hard if you try (except maybe the Kaiser Window stuff). Anyway, here is a simple JS implementing the filter:

desc:Linear-phase FIR filter
// (c) Theo Niessink 2013
// License: GPL - http://www.gnu.org/licenses/gpl.html

slider1:1<0,4,1{Bypass,Low-Pass,High-Pass,Band-Pass,Band-Stop}>Filter Type
slider2:1000<20,20000,1>Frequency (Hz)
slider3:2.0<0.0,4.0,0.1>Bandwidth (oct)
slider4:63<1,512,1>Length (samples)
slider5:3<0,5,1{Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser}>Window Type

// Source: http://www.labbookpages.co.uk/audio/firWindowing.html

@init

pdc_bot_ch = 0;
pdc_top_ch = 2;

function fact(x)
  local(y)
(
  y = (x |= 0);
  while(x > 1 ? y *= (x -= 1));
  max(y, 1);
);

function i0(x)
  local(y, i)
(
  y = i = 1;
  loop(20,
    y += (x/2)^(2*i) / fact(i)^2;
    i += 1;
  );
  y;
);

alpha = 3.0;

@slider

filter = slider1|0;
window = slider5|0;

ft = min(slider2 / srate, 0.5);
filter >= 3 ? (
  ft1 = ft * 2^(-0.5 * slider3);
  ft2 = ft * 2^( 0.5 * slider3);
);

len = slider4|0;
filter > 1 ? len |= 1;
m = len - 1;

n = 0;
loop(len,

  filter == 1 ? w[n] = n != m/2 ?  sin(2*$pi * ft  * (n - m/2)) / ($pi * (n - m/2))  :     2 *  ft :
  filter == 2 ? w[n] = n != m/2 ? -sin(2*$pi * ft  * (n - m/2)) / ($pi * (n - m/2))  : 1 - 2 *  ft :
  filter == 3 ? w[n] = n != m/2 ?  sin(2*$pi * ft2 * (n - m/2)) / ($pi * (n - m/2)) -
                                   sin(2*$pi * ft1 * (n - m/2)) / ($pi * (n - m/2))  :     2 * (ft2 - ft1) :
  filter == 4 ? w[n] = n != m/2 ?  sin(2*$pi * ft1 * (n - m/2)) / ($pi * (n - m/2)) -
                                   sin(2*$pi * ft2 * (n - m/2)) / ($pi * (n - m/2))  : 1 - 2 * (ft2 - ft1);

  window == 0 ? w[n] *= 1 :
  window == 1 ? w[n] *= 1 - 2 * abs(n - m/2) / m :
  window == 2 ? w[n] *= 0.5 - 0.5 * cos(2*$pi*n / m) :
  window == 3 ? w[n] *= 0.54 - 0.46 * cos(2*$pi*n / m) :
  window == 4 ? w[n] *= 0.42 - 0.5 * cos(2*$pi*n / m) + 0.08 * cos(4*$pi*n / m) :
  window == 5 ? w[n] *= i0($pi * alpha * sqrt(1 - (2*n / m - 1)^2)) / i0($pi * alpha);

  n += 1;
);

buf = len;
pdc_delay = len / 2 |0;

@block

s = buf + samplesblock;
memcpy(s, buf, m);

@sample

s -= 1;
s[] = 0.5 * (spl0 + spl1);

filter ? (
  sum = n = 0;
  loop(len,
    sum += w[n] * s[n];
    n += 1;
  );
  spl0 = sum;
) : (
  spl0 = s[pdc_delay];
);
spl1 = spl0;

[IXix]

Originally Posted by Tale

Well, actually it isn't that hard if you try

I think you overestimate my intelligence!

Seriously though, I think I'm almost starting to understand how it works. Example code is always helpful so thanks again.

[Tale]

Originally Posted by Tale

Well, actually it isn't that hard if you try (except maybe the Kaiser Window stuff).

I have tried a little harder myself, and I have now added the Kaiser window.

EDIT: And I have slightly improved CPU use by moving memcpy() to @block and adding s[].

[Tale]

Here is a revamped, more practical windowed FIR filter example. Changes:

• Predefined filter length, which is scaled to sample rate, so the filter sounds the same at different sample rates.
• Filter length is always odd.
• Predefined Blackman window, which is precalculated in @init.
• Because both the filter and window are symmetrical, only half is actually caclulated.
• Optimized LPF, HPF, BP, BS, and window caclulations.

desc:Linear-phase FIR filter
// (c) Theo Niessink 2013
// License: GPL - http://www.gnu.org/licenses/gpl.html

slider1:1<0,4,1{Bypass,Low-Pass,High-Pass,Band-Pass,Band-Stop}>Filter Type
slider2:1000<20,20000,1>Frequency (Hz)
slider3:2.0<0.0,4.0,0.1>Bandwidth (oct)

// Source: http://www.labbookpages.co.uk/audio/firWindowing.html

@init

len = 63 / 44100 * srate | 1;
m = len - 1;

pdc_bot_ch = 0;
pdc_top_ch = 2;
pdc_delay = len / 2 | 0;

wnd = 0;
w = wnd + pdc_delay;
buf = w + len;

function blackman(n, m)
  local(x)
(
  x = 2*$pi*n / m;
  0.42 - 0.5 * cos(x) + 0.08 * cos(2*x);
);

n = 0;
loop(pdc_delay,
  wnd[n] = blackman(n, m);
  n += 1;
);

function lpf(n, m, ft)
  local(x)
(
  x = $pi * (n - m/2);
  sin(2*ft*x) / x;
);

function hpf(n, m, ft)
(
  -lpf(n, m, ft);
);

function bp(n, m, ft1, ft2)
  local(x)
(
  x = $pi * (n - m/2);
  (sin(2*ft2*x) - sin(2*ft1*x)) / x;
);

function bs(n, m, ft1, ft2)
  local(x)
(
  x = $pi * (n - m/2);
  (sin(2*ft1*x) - sin(2*ft2*x)) / x;
);

@slider

filter = max(0, slider1 + 0.5 |0);
bw = 2^(0.5 * slider3);

ft = min(0.5, slider2 / srate);
ft1 = min(0.5, ft / bw);
ft2 = min(0.5, ft * bw);

filter == 1 ? (
  n = 0;
  loop(pdc_delay,
    w[n] = wnd[n] * lpf(n, m, ft);
    n += 1;
  );
  w[n] = 2*ft;
) :

filter == 2 ? (
  n = 0;
  loop(pdc_delay,
    w[n] = wnd[n] * hpf(n, m, ft);
    n += 1;
  );
  w[n] = 1 - 2*ft;
) :

filter == 3 ? (
  n = 0;
  loop(pdc_delay,
    w[n] = wnd[n] * bp(n, m, ft1, ft2);
    n += 1;
  );
  w[n] = 2*(ft2 - ft1);
) :

filter == 4 ? (
  n = 0;
  loop(pdc_delay,
    w[n] = wnd[n] * bs(n, m, ft1, ft2);
    n += 1;
  );
  w[n] = 1 - 2*(ft2 - ft1);
);

filter ? loop(pdc_delay,
  n += 1;
  w[n] = w[m - n];
);

@block

s = buf + samplesblock;
memcpy(s, buf, m);

@sample

s -= 1;
s[] = 0.5 * (spl0 + spl1);

filter ? (
  sum = n = 0;
  loop(len,
    sum += w[n] * s[n];
    n += 1;
  );
  spl0 = sum;
) : (
  spl0 = s[pdc_delay];
);
spl1 = spl0;

[IXix]

Thanks Tale! Much appreciated!

[IXix]

Originally Posted by Tale

pdc_delay = len / 2 | 0;

Why is pdc_delay only half the filter length? I can see that value is used in calculating the window/filter weights because you only need half the values but I don't understand why pdc_delay isn't the full filter length.

I think I'm nearly, almost, not quite, just about starting to get this.

[bang]

Originally Posted by IXix

Why is pdc_delay only half the filter length?

because the output of the filter is at the midpoint of the sync window. the filter is calculated using samples before and after the output sample. but the pdc_delay is only the after samples.

sorry to have been awol on this. life's been a bit overwhelming. big thanks to Tale for jumping in with code and the excellent labbookpages link. i will still do some comments/sample code for the optimizations rsn, since that was a pretty big win.

enjoy! /dan

[IXix]

Originally Posted by bang

because the output of the filter is at the midpoint of the sync window. the filter is calculated using samples before and after the output sample. but the pdc_delay is only the after samples.

Ahhh! Got it.

Originally Posted by bang

sorry to have been awol on this. life's been a bit overwhelming. big thanks to Tale for jumping in with code and the excellent labbookpages link. i will still do some comments/sample code for the optimizations rsn, since that was a pretty big win.

I'd still like a simple example of how to use your include but no worries, life can really get in the way of being a nerd sometimes. In the meantime, maybe I'll be able to figure it out myself now I'm starting to see how the whole thing works.

[mwe]

After days of banging my head against this I think I finally understand how it works. Thanks for the references and sample code Tale, I always feel like I take a major step forward when I dissect your stuff.

[IXix]

I think my brain is having an off day. How would I apply this to a chunk of audio in memory rather than a stream of incoming samples? I *think* I just have to move the window position along one sample at a time and overwrite the samples as they're processed.

If I'm right about that, how would the audio start/end points be handled?

At the start of the audio, would the window be aligned with the first sample or would it need to precede the audio by filter length - 1?

Start here, window aligned with the first sample?
[.........]  <- Stored samples
[...]        <- Filter window

or here, window overlapping only the first sample?
  [.........]  <- Stored samples
[...]          <- Filter window

And at the end I guess I stop when the filter window is aligned with the last sample but that means that the window will be mostly over nothing. Will that work or would I have to stop when the window goes over the end of the audio (ie. don't process silence)

Stop here, while there are still samples left to process?
[.........]  <- Stored samples
      [...]  <- Filter window

or here, allowing the window to process some nothingness?
[.........]   <- Stored samples
        [...] <- Filter window

[Tale]

I don't think you necessarily have to to anything special at the beginning. At the end I guess you will have to continue until you have outputed the number of input samples plus the length of the buffer.

But of course then you will end up with an audio signal that is delayed by half the filter length, so... I guess you might want to cut off this amount of samples from the beginning after all. In that case you would end up with the number of input samples plus half the filter length.

[IXix]

Originally Posted by Tale

I don't think you necessarily have to to anything special at the beginning. At the end I guess you will have to continue until you have outputed the number of input samples plus the length of the buffer.

But of course then you will end up with an audio signal that is delayed by half the filter length, so... I guess you might want to cut off this amount of samples from the beginning after all. In that case you would end up with the number of input samples plus half the filter length.

Thanks Tale. The audio isn't actually destined to be output so there's no problem with delay, although I may need to discard the extra samples. I just need to do a lowpass on the samples before analysing them with some other code but since I won't be able to hear it, I'm not sure how to check lowpass results, which is why I need to be sure that the method is correct. I can't just hack at it until it sounds right.

Working on something cool for an extension.

[SaulT]

I really like the code, it is more compact and easier to understand (and includes highpass/bandpass) than what I've written. I'm not seeing coefficient normalization though. I don't know if that's something that people necessarily are looking for, but it was crucial for me when doing oversampling.

The sum of the coefficients of a FIR is the multiplier of the signal level (DC gain?). So, if they all add to 1, then the signal level would be unchanged. Here's one way to do it - first accumulating the coefficient values, then multiplying them by the normalizing factor. Because you're only calculating half of the values, all but the middle value should be doubled.

n = 1;
normalize = wdn[0];
loop(pdc_delay-1,
  normalize += wdn[n] * 2; // we have to account for the mirrored side, too
  n += 1;
);

n = 0;
normalize = 1/normalize;
loop(pdc_delay,
  wdn[n] *= normalize;
  n += 1;
);

The problem that I ran into was that at srate/2 the audio level wasn't changed, but as I lowered the ft it lost a lot of signal... which makes sense, but it's not necessarily ideal, especially if you're doing something like oversampling.

This is a link to my windowed sinc lowpass plugin. The code isn't as pretty or as functional as Tale's in some ways, although FWIW I do have some window methods that he doesn't.

http://forum.cockos.com/showthread.php?t=136765

[SaulT]

or here, allowing the window to process some nothingness?
[.........]   <- Stored samples
        [...] <- Filter window

The purpose of the pdc_delay is to account for this - it delays the results by the proper number of samples for you. If you were not to use pdc_delay, that's when you'd start having to account for the delays. If you are processing a continuous stream of audio, using pdc_delay is the proper method.

*** technical stuff

Window methods and tap counts are trade-offs. The Blackman window has the least spectral leakage of the window methods that I've tried, but it's transient response is by comparison so-so. The Flat-top method has much better transient response, but more spectral leakage. The more taps you use the better the resolution (the narrower the passband), but the more likely you are to have pre-ringing. It's a trade-off between performance in the frequency domain and frequency in the time domain. Like, there are no free lunches.

In other words, a tap count in the thousands will give you a few milliseconds of pre-ringing, while a tap count in the tens will give you practically none. A tap count in the thousands will have a bandwidth (width of passband) in the tens of hertz, while a tap count in the tens will have a passband width in the thousands of hertz.

Maybe I can give an example to illustrate? Let's imagine a FIR filter with the coefficients [0.1 0.2 0.4 0.2 0.1]... kinda like a triangular window:

0.1  0.2  0.4  0.2  0.1   output

0    0    0    0    0     0
0    0    0    0    1     0.1  <-- pre-ringing
0    0    0    1    0     0.2
0    0    1    0    0     0.4  <-- sample delay 
0    1    0    0    0     0.2  <-- post-ringing
1    0    0    0    0     0.1  
0    0    0    0    0     0

We see that as soon as a sample enters the FIR the FIR begins to output something (usually a small value). After the delay it still outputs a small value. This is a slightly exaggerated example, in real life a windowed sinc filter will have much small values at the edges and so have very little pre- and post-ringing (by comparison), but it will still be there.

...

Umm, that's kind've a tangent. Maybe it's useful?

[IXix]

All information is useful, so thanks for contributing to my limited understanding of filters. I've had to put my idea on hold but I'll get back to it eventually.

[Smashed Transistors]

If you want to complement FIRs with some IIR filters, you can have a look to Bessel filters. These iir filters have a maximally flat group delay response. Can be used along with delays to have complementary HP LP filters for frequency splitting.

[Tale]

Originally Posted by Smashed Transistors

If you want to complement FIRs with some IIR filters, you can have a look to Bessel filters. These iir filters have a maximally flat group delay response. Can be used along with delays to have complementary HP LP filters for frequency splitting.

Yeah, Bessel IIR filters are nice because they are much cheaper to process (when compared to brute force FIR), especially at high oversampling factors, because with FIR you need to scale the kernel length by the oversampling factor, so brute force quickly becomes impractical. However, caclulating Bessel IIR filter coefficients is hard work (that is, for the programmer, not so much for the CPU). Anyway, there is a Bessel filter implementation in my WDL.git repository, but it's in C++, not (yet) in JS...

--

Anyway, I have been looking into the overlap-add method to apply convolution using FFT, so I thought I could reimplement my FIR filter example using that:

desc:Linear-phase FIR filter
// (c) Theo Niessink 2013, 2014
// License: GPL - http://www.gnu.org/licenses/gpl.html

slider1:1<0,4,1{Bypass,Low-Pass,High-Pass,Band-Pass,Band-Stop}>Filter Type
slider2:1000<20,20000,1>Frequency (Hz)
slider3:2.0<0.0,4.0,0.1>Bandwidth (oct)

// Sources:
// http://www.labbookpages.co.uk/audio/firWindowing.html
// http://en.wikipedia.org/wiki/Overlap-add_method

@init

function int(x) ( x|0 );
function log2(n) ( ceil(log(n) / log(2)) );

// FIR length needs to be odd for HP/BP/BS.
fir_len = int(63 / 44100 * srate) | 1;

fft_len = 1 << log2(2*fir_len - 1);
fft_len > 32768 ? (
  fft_len = 32768;
  fir_len = fft_len/2 - 1;
);
scale = 1/fft_len;

m = fir_len - 1;

pdc_bot_ch = 0;
pdc_top_ch = 2;
pdc_delay = int(1.5*fir_len);

h = 0;
in = h + 2*fft_len;
out = in + 2*fft_len;

function blackman(n, m)
  local(x)
(
  x = 2*$pi*n / m;
  0.42 - 0.5 * cos(x) + 0.08 * cos(2*x);
);

function lpf(n, m, ft)
  local(x)
(
  x = $pi * (n - m/2);
  sin(2*ft*x) / x;
);

function hpf(n, m, ft)
(
  -lpf(n, m, ft);
);

function bp(n, m, ft1, ft2)
  local(x)
(
  x = $pi * (n - m/2);
  (sin(2*ft2*x) - sin(2*ft1*x)) / x;
);

function bs(n, m, ft1, ft2)
  local(x)
(
  x = $pi * (n - m/2);
  (sin(2*ft1*x) - sin(2*ft2*x)) / x;
);

@slider

filter = max(0, int(slider1 + 0.5));
bw = 2^(0.5 * slider3);

ft = min(0.5, slider2 / srate);
ft1 = min(0.5, ft / bw);
ft2 = min(0.5, ft * bw);

n = 0;
half = int((m + 1) / 2);

filter == 1 ? (
  loop(half,
    h[2*n] = blackman(n, m) * lpf(n, m, ft);
    h[2*n + 1] = 0;
    n += 1;
  );
  h[2*n] = 2*ft;
) :

filter == 2 ? (
  loop(half,
    h[2*n] = blackman(n, m) * hpf(n, m, ft);
    h[2*n + 1] = 0;
    n += 1;
  );
  h[2*n] = 1 - 2*ft;
) :

filter == 3 ? (
  loop(half,
    h[2*n] = blackman(n, m) * bp(n, m, ft1, ft2);
    h[2*n + 1] = 0;
    n += 1;
  );
  h[2*n] = 2*(ft2 - ft1);
) :

filter == 4 ? (
  loop(half,
    h[2*n] = blackman(n, m) * bs(n, m, ft1, ft2);
    h[2*n + 1] = 0;
    n += 1;
  );
  h[2*n] = 1 - 2*(ft2 - ft1);
);

filter ? (
  !(m & 1) ? (
    h[2*n + 1] = 0;
    n += 1;
  );
  loop(half,
    h[2*n] = h[2*(m - n)];
    h[2*n + 1] = 0;
    n += 1;
  );
);

memset(h + 2*fir_len, 0, 2*(fft_len - fir_len));
fft(h, fft_len);

@sample

in[2*idx] = 0.5 * (spl0 + spl1);
in[2*idx + 1] = 0;
idx += 1;

filter ? (
  idx >= fir_len ? (
    memset(in + 2*idx, 0, 2*(fft_len - idx));
    idx = 0;

    fft(in, fft_len);
    // WTF, there's complex multiplication built into JSFX?! Cool!!1!
    convolve_c(in, h, fft_len);
    ifft(in, fft_len);

    memcpy(out, out + fir_len, fft_len - fir_len);
    memset(out + fir_len, 0, fft_len - fir_len);

    n = 0;
    loop(2*fir_len - 1,
      out[n] += in[2*n] * scale;
      n += 1;
    );
  );

) :

/* !filter ? */ (
  idx >= pdc_delay ? (
    idx = 0;

    n = 0;
    loop(pdc_delay,
      out[n] = in[2*n];
      n += 1;
    );
  );
);

spl0 = spl1 = out[idx];

This is faster than my earlier brute force method (about 4x on my Core i7), even at 64 samples or less. And the beauty is that if you increase the FIR kernel length, it doesn't really become much slower. Now, the brute force method will be almost N times slower when the FIR kernel is N times longer.

There is a catch however (isn't there always): The latency of the brute force method is "only" half the FIR kernel, but using this method the latency is increased with half the FFT size, so this means the latency will be about 1.5x the FIR kernel length.

Note that for simplicity I do only a single fft > convolve_c > ifft per buffered sample block, so the maximum FIR kernel length (16384) is half the maximum FFT (32768). I guess next up would be to split up longer kernels into smaller FFT chunks, which should remove this limit.

EDIT: Fixed FIR length, needs to be odd for HP/BP/BS.
EDIT: Minor fixes (buffer fir_len not fft_len/2, latency is now exactly 1.5*fir_len, fixed bypass).

[bertrand fraysse]

wow this is really great ! Thank you.
I spotted a mistake on the latency calculation, it has to be

pdc_delay = int(1.5*fir_len)-1;

to perfectly null the output when in delta solo mode, and produce phase consistency with unfiltered signals.

I'll try and make everything in functions to be able to call multiple filters at the same time.

Thanks again for the good and clear work !

Maybe you expand on this, as it's been almost 10 years ago ?

EDIT = I'm having difficulties to understand the memset, memcopy and what to do if I want a stereo version of these filters
I'd like to try to make a graphic EQ out of it too but it's out of my knowledge. if anyone could shime in it would be great ! thanks.

[bertrand fraysse]

I managed to encapsulate everything into functions
and made it work.
I had problems with memory managment.

I tried to simplify the code, and hopefully minimize the CPU use.
I hope you'll enjoy.

Would there be a way to avoid the clicks when moving the frequency cutoff ?
I'm pretty sure, no as interpolation of all the calculations would be crazy amount of CPU.

desc:Linear-phase FIR filter
// (c) Theo Niessink 2013, 2014 with some modifications by Bertrand Fraysse
// License: GPL - http://www.gnu.org/licenses/gpl.html

slider1:1<0,4,1{Bypass,Low-Pass,High-Pass,Band-Pass,Band-Stop}>Filter Type
slider2:0<0,1,.0001>Frequency (Hz)
slider3:1<0.0,16,0.1>Bandwidth (oct)

// Sources:
// http://www.labbookpages.co.uk/audio/firWindowing.html
// http://en.wikipedia.org/wiki/Overlap-add_method

@init

twopi = $pi+$pi;
function int(x) ( x|0 );
function log2(n) ( ceil(log(n) / log(2)) );

// FIR length needs to be odd for HP/BP/BS.
fir_len_ = 61;

pdc_bot_ch = 0;
pdc_top_ch = 2;
pdc_delay = int(1.5*fir_len_)-1;

function filterInit (fir_len_, memIncrement_)
instance (memIncrement, fir_len, two_fir_len, fft_len, scale, two_fft_len, fft_len_minus_fir_length,
two_fft_len_minus_fir_length, m, half_m, h, in, out)
(
  fir_len = fir_len_;
  memIncrement = memIncrement_;
  fft_len = 1 << log2(2*fir_len - 1);
  fft_len > 32768 ? (
    fft_len = 32768;
    fir_len = fft_len/2 - 1;
  );
  
  two_fir_len = fir_len+fir_len;
  scale = 1/fft_len;
  two_fft_len = fft_len+fft_len;
  fft_len_minus_fir_length = fft_len - fir_len;
  two_fft_len_minus_fir_length = fft_len_minus_fir_length+fft_len_minus_fir_length;

  m = fir_len - 1;
  half_m = m*.5;

  h = 0+memIncrement;
  in = h + two_fft_len;
  out = in + two_fft_len;
);

function blackman(n)
  local(x)
  instance (m)
  global (twopi)
(
  x = twopi*n / m;
  .42 - .5 * cos(x) + .08 * cos(x+x);
);

function lpf(n, ft)
  local(x)
  instance (half_m)
(
  x = $pi * (n - half_m);
  sin(ft*x) / x;
);

function hpf(n, ft)
(
  -this.lpf(n, ft);
);

function bp(n, ft1, ft2)
  local(x)
  instance (half_m)
(
  x = $pi * (n - half_m);
  (sin(ft2*x) - sin(ft1*x)) / x;
);

function bs(n, ft1, ft2)
  local(x)
  instance (half_m)
(
  x = $pi * (n - half_m);
  (sin(ft1*x) - sin(ft2*x)) / x;
);

function filterCoefCalculation (filter, ft, bw)
instance (two_ft, ft1, ft2, two_ft1, two_ft2, two_bandsize, n, m, half, half_m, h, two_fft_len_minus_fir_length, two_fir_len, fft_len)
(
  ft = ft*ft*ft*ft*.5;
  two_ft = ft+ft;
  ft1 = min(0.5, ft / bw);
  ft2 = min(0.5, ft * bw);
  two_ft1 = ft1+ft1;
  two_ft2 = ft2+ft2;
  two_bandsize = ft2-ft1 + ft2-ft1;
  
  n = 0;
  half = int((m + 1)*.5);
  
  filter == 1 ? (
    loop(half,
      h[n+n] = this.blackman(n) * this.lpf(n, two_ft);
      h[n+n + 1] = 0;
      n += 1;
    );
    h[n+n] = ft+ft;
  ) :
  
  filter == 2 ? (
    loop(half,
      h[n+n] = this.blackman(n) * this.hpf(n, two_ft);
      h[n+n + 1] = 0;
      n += 1;
    );
    h[n+n] = 1 - 2*ft;
  ) :
  
  filter == 3 ? (
    loop(half,
      h[n+n] = this.blackman(n) * this.bp(n, two_ft1, two_ft2);
      h[n+n + 1] = 0;
      n += 1;
    );
    h[n+n] = two_bandsize;
  ) :
  
  filter == 4 ? (
    loop(half,
      h[n+n] = this.blackman(n) * this.bs(n, two_ft1, two_ft2);
      h[n+n + 1] = 0;
      n += 1;
    );
    h[n+n] = 1 - two_bandsize;
  );
  
  filter ? (
    !(m & 1) ? (
      h[n+n + 1] = 0;
      n += 1;
    );
    loop(half,
      h[n+n] = h[2*(m - n)];
      h[n+n + 1] = 0;
      n += 1;
    );
  );
  memset(h + two_fir_len, 0, two_fft_len_minus_fir_length);
  fft(h, fft_len);
);

function filterProcess (input)
instance (memIncrement, in, output, idx, h, fir_len, two_fft_len_minus_idx, fft_len, out, fir_len, fft_len_minus_fir_length, two_fir_len, n2, scale)
(
in[idx+idx] = input;
in[idx+idx + 1] = 0;

idx += 1;
  idx >= fir_len ?
  (
    two_fft_len_minus_idx = 2*(fft_len - idx);
    memset(in + idx+idx, 0, two_fft_len_minus_idx);
    idx = 0;

    fft(in, fft_len);
    convolve_c(in, h, fft_len);
    ifft(in, fft_len);

    memcpy(out, out + fir_len, fft_len - fir_len);
    memset(out + fir_len, 0, fft_len_minus_fir_length);

    n2 = 0;
    loop(two_fir_len - 1,
      out[n2] += in[n2+n2] * scale;
      n2 += 1;
    );
  );
  output = out[idx];
);

filter1.filterInit (fir_len_, 0);
filter2.filterInit (fir_len_, filter1.out*2);

@slider

  filter = max(0, int(slider1 + 0.5));
  bw = 2^(0.5 * slider3);
  ft = slider2;
  
filter1.filterCoefCalculation (filter, ft, bw);
filter2.filterCoefCalculation (filter, ft, bw);


@sample

input1 = spl0;
input2 = spl1;

filter1.filterProcess (input1);
filter2.filterProcess (input2);

spl0 = filter1.output;
spl1 = filter2.output;

[Tale]

Originally Posted by bertrand fraysse

I spotted a mistake on the latency calculation, it has to be

pdc_delay = int(1.5*fir_len)-1;

to perfectly null the output when in delta solo mode, and produce phase consistency with unfiltered signals.

You're right! However, I think the FFT introduces fir_len-1 samples of latency, and the FIR filter itself (fir_len-1)/2, so I guess it should be:

pdc_delay = int(1.5*(fir_len - 1));

But I haven't checked if this is actually correct.

Originally Posted by bertrand fraysse

Would there be a way to avoid the clicks when moving the frequency cutoff ?

Maybe you could run two filters in parallel, and then crossfade between them?

[bertrand fraysse]

Both modifications of the pdc_delay exhibit the same behaviors, so I'll use yours.

Originally Posted by Tale

Maybe you could run two filters in parallel, and then crossfade between them?

I just thought about it this morning :-)
I've made some crossfading delays that sound great, I may adapt the same algorithm at some point, in this case it's morce complicated with the coef calculation, but possible anyway.

Thank you again ! I really like those filters, they are pretty cheap on the CPU, the slope can be very steep, and they introduce no phase issues. what may be the downside of those ? maybe not the best in the low frequencies...

[bertrand fraysse]

I tried with my crossfading algorithm but my concept is not working here
as the coeff calculation takes time, I get a click when my crossfading is done.
my algorithm is

if there's new coeff
create a new filter with the new coef and crossfade it with the previous filter.
when the crossfading is done, actualise the previous filter coeff.

maybe I should make some sort of pingpong between the two filters.
I may try that, but if you already have a nice way to do this, I may be interested in having a look.
thank you.

[Mikobuntu]

I tested @bertrand fraysse plugin. I ran it through Bertom EQ curve Analyzer and noticed that the filters for example using high pass there is a low shelf cut applied up to around 1khz. Was just wondering is this the expected behaviour ?
thanks guys. [[ see attached .gif ]] Sorry for the bad quality the upload size limit required me to scale this down a lot.

[bertrand fraysse]

Originally Posted by Mikobuntu

I tested @bertrand fraysse plugin. I ran it through Bertom EQ curve Analyzer and noticed that the filters for example using high pass there is a low shelf cut applied up to around 1khz. Was just wondering is this the expected behaviour ?

you may want to try increasing the parameter -> fir_len_ = 61; you can go pretty high without adding too much CPU, something like 255 ?

Anyway, this may explain what I felt on the low frequency range.
Seems to be an interesting curve...

did not knew about this Bertom EQ curve analyzer. this seem really useful.
Thank you.

[KostyaK]

Hello
can your linear-phase filter be used to make a mono-maker?