Linear phase filter

trogluddite · 9 Jul 2012

Warmonger wrote: the design algorithms I found need quite complex complex (sic!) math and multiple iterations

Unfortunately, that is because those algorithm are necessary to make linear phase work. Almost all SM filters are based around the bi-quad algorithm, because it is very light on CPU. That design is known as an Infinite Impulse Response filter (IIR) - Infinite, because it is effectively two very short delay lines with feedback; so there's always a tiny amount of ringing that never quite dies away (or worse, turns into the dreaded de-normals).

To make a linear phase filter you have to use the alternative Finite Impulse Response algorithm (FIR) - which is more like a delay with LOTS of taps, but no feedback (technically it's called convolution - exactly the same as one of those impulse Response type reverb plugins). In principle, simple enough, but updating all those delay taps is very CPU heavy, so many designs use a faster method related to FFT's - which is where all the gobbledy-gook complex maths comes from.

To get an FIR filter to be linear phase, there's also another condition - half of the taps need to be PRE-delays that hear the audio before the current sample position. In fact what really happens is that the audio is delayed slightly; the host can compensate for this latency on Playback (but not live Recording) by using Plugin Delay Compensation.

If you have your heart set on perfect linear phase, search the forum etc. for some of those techy words like 'convolution' - I'm pretty sure that there have been some experiments with it before. The Wikipedia pages aren't too bad either, if you want to understand more about IIR/FIR differences. Whether you can get even close with regular IIR filters, I'm not sure (hopefully, Martin will be around shortly!) - no doubt some must be better than others, or maybe there is some trickery with all-pass filters that can be done to straighten out the phase.

martinvicanek · 10 Jul 2012

Warmonger, what other characteristics besides linear phase should the filter have? What specific application do you have in mind?

it can be just one-pole

That's a bit of a contradiction: poles occur in IIR (= recursive) filters, which do not have linear phase.

Warmonger · 10 Jul 2012

I want to put it in a feedback loop to create 'airy' reverb. Linear phase is necessary to keep constant echo density of reverb, but attenuation can be minimal.

oddson · 10 Jul 2012

I once experimented with FIR filters.

If symetricallity of coefficients is a sufficient condition for linear phase response (I think it is but maybe Trog or another DSP person can confirm) then this file might be helpful in playing with coefficients.

The 'Smoother Control' works as a xy pad to control the coefficents --- just play with it, its simpler than explaining.

(The weighted moving average module takes a scalar value and applies an exponential weighting to generate the coeffeicients -- this is certainly not an 'optimal' method but it's what I was able to generate on my own with my limited understanding and it seems to work and it doesn't require the complex complex maths that are required when optimality is sought.)

This was a VERY long time ago so I'm not sure about the guts of it any more but maybe you'll find it helpful.

oddson · 11 Jul 2012

I can't find a 'definitive' source right now but a couple of blog posts from smarty-pants DSP types suggest coefficients that are symmetrical around the centre is both necessary and sufficient for linear phase.

If you don't want to adjust the response curve dynamically you can use some of the fancy optimal calculation techniques like "Remez exchange" with an online calculator that spit out coefficients and just plug them in manually. You don't need to understand the math.

cyto · 11 Jul 2012

Warmonger wrote: I want to put it in a feedback loop to create 'airy' reverb. Linear phase is necessary to keep constant echo density of reverb, but attenuation can be minimal.

Hmmmm..... I'm not sure that going for a completely phase-accurate response is necessarily a good idea in this situation. I'd be very curious to see what algorithm you are using to create the reverb. Static echo density will no doubt provide a resonant mess, with that "metallic" sound that plagues a lot of digital reverbs. Reverb has been somewhat of an obsession of mine in recent weeks and I have been in the process of recreating some of the more common algorithms floating around in an attempt to find a good solution to replace the stock SM reverb (which I hate!). Most of the more successful algorithms (lexicon, eventide, etc.) aim for the echo density to "build" (like it would in a real room). Even the "old standard" Schroeder algorithms strive to "smear" the phase response (In fact, they are essentially just a network of combs and APs).

Having said that, don't discount just the standard first-order HP IIR filter. As you can see, the phase response is fairly flat within only a couple octaves over the corner frequency...

Notice that this is tuned to 160Hz. Within the first octave, it is only about 10-15 degrees out of phase. (the white line is the phase offset and is centered around zero, the red line is the magnitude with 0dBFS up at the top). If you want to try this filter and see how it sounds, the transfer function is...

streamin in; streamin hpCoeff; streamout out;

float in1; out = in - in1 + (hpCoeff * out); in1 = in;

You can calculate the coefficient in green as hpCoeff=e^(-pi *cornerFreq/sampleFreq). where e = 2.7182817 and -pi = -3.1415927

I imagine that you have many delay lines going and to run each one through a FIR is going to kill your CPU. I worry that all that processing might negatively affect your reverb by creating more resonant ringing. Try the first order and see if it's okay (it's very light). But again, if your algorithm definitely calls for perfect phase accuracy, then by all means go for it! I don't want to discourage you! In fact, I'd love to look at the algorithm myself. 🙂

-cyto

Warmonger · 11 Jul 2012

In fact, I tried to using feedback delay network. After making sure it sounds perfect (no metallic ringing) I realized it's not exactly FDN according to schematic, but still works nicely.

EDIT: It works with your filter, thanks. The air effect itself may be not than noticeable, but certainfly filtering lower frequencies dramatically improved the quality of sound. So natural 🙂

infuzion · 25 Jul 2012

I also am looking for a Linear phase filter, though for a frequency splitter. I think cyto's example will cause too much phasing when stacked to get the cut angle I want for a splitter. Ideas please?

martinvicanek · 28 Jul 2012

Here is a basic windowed sinc FIR filter. It is linear phase since the filter coefficients are symmetric wrt the center. It will do OK for cutoffs in the kHz range - for lower cutoffs you will need to increase the filter order (31 in this example).

Here is a good overview on FIR filter design.

infuzion · 28 Jul 2012

martinvicanek wrote: Here is a basic windowed sinc FIR filter. It is linear phase since the filter coefficients are symmetric wrt the center. It will do OK for cutoffs in the kHz range - for lower cutoffs you will need to increase the filter order (31 in this example).

Here is a good overview on FIR filter design.

Both links are great; thanks! Sorry, I was not specific about my requirements; I need to low-pass base signals, likely even to 15 or 20Hz (mastering). Do I add more "n" in your formula? TIA

BTW, I would pre-calculate n - 16; though not much of a speed increase on its own due to the hop(), it has additional benefits, like buffering, helping with ASM compiler optimization, etc. n = 1; nM16 = n-16; w = 1 + cos1(nM16/32); c1 = w*sin1(nM16*fh)/nM16; c31 = c1 ; n = 2; nM16 = n-16; w = 1 + cos1(nM16/32); c2 = w*sin1(nM16*fh)/nM16; c30 = c2 ; n = 3; nM16 = n-16; w = 1 + cos1(nM16/32); c3 = w*sin1(nM16*fh)/nM16; c29 = c3 ; n = 4; nM16 = n-16; w = 1 + cos1(nM16/32); c4 = w*sin1(nM16*fh)/nM16; c28 = c4 ; n = 5; nM16 = n-16; w = 1 + cos1(nM16/32); c5 = w*sin1(nM16*fh)/nM16; c27 = c5 ; n = 6; nM16 = n-16; w = 1 + cos1(nM16/32); c6 = w*sin1(nM16*fh)/nM16; c26 = c6 ; n = 7; nM16 = n-16; w = 1 + cos1(nM16/32); c7 = w*sin1(nM16*fh)/nM16; c25 = c7 ; n = 8; nM16 = n-16; w = 1 + cos1(nM16/32); c8 = w*sin1(nM16*fh)/nM16; c24 = c8 ; n = 9; nM16 = n-16; w = 1 + cos1(nM16/32); c9 = w*sin1(nM16*fh)/nM16; c23 = c9 ; n = 10; nM16 = n-16; w = 1 + cos1(nM16/32); c10 = w*sin1(nM16*fh)/nM16; c22 = c10; n = 11; nM16 = n-16; w = 1 + cos1(nM16/32); c11 = w*sin1(nM16*fh)/nM16; c21 = c11; n = 12; nM16 = n-16; w = 1 + cos1(nM16/32); c12 = w*sin1(nM16*fh)/nM16; c20 = c12; n = 13; nM16 = n-16; w = 1 + cos1(nM16/32); c13 = w*sin1(nM16*fh)/nM16; c19 = c13; n = 14; nM16 = n-16; w = 1 + cos1(nM16/32); c14 = w*sin1(nM16*fh)/nM16; c18 = c14; n = 15; nM16 = n-16; w = 1 + cos1(nM16/32); c15 = w*sin1(nM16*fh)/nM16; c17 = c15;

Warmonger · 28 Jul 2012

Nice try, but...

martinvicanek · 28 Jul 2012

infuzion wrote: I need to low-pass base signals, likely even to 15 or 20Hz (mastering). Do I add more "n" in your formula?

Wow, that's deep, not even Trog's low B would pass such a LP, hehe 😉 . For such a low cutoff a FIR filter is very inefficient, you'd need thousands of operations per sample. In that case perhaps an IIR filter with subsequent allpass phase equalization in the passband as Trog suggested would be better.

infuzion wrote: BTW, I would pre-calculate n - 16; though not much of a speed increase on its own due to the hop(), it has additional benefits, like buffering, helping with ASM compiler optimization, etc.

Oh, sure, there is plenty of room for optimization! For instance it is a complete waste to use sin1() and cos1(), much better to compute recursively.

infuzion · 28 Jul 2012

martinvicanek wrote: infuzion wrote: I need to low-pass base signals, likely even to 15 or 20Hz (mastering). Do I add more "n" in your formula?

Wow, that's deep, not even Trog's low B would pass such a LP, hehe 😉 . For such a low cutoff a FIR filter is very inefficient, you'd need thousands of operations per sample. In that case perhaps an IIR filter with subsequent allpass phase equalization in the passband as Trog suggested would be better.

Well, for mastering 1000s of ops would be OK. But for my immediate application I think down to 100Hz or so would be fine.

Seems that warmonger found a better phase-finding tool... maybe another try please? 😇

martinvicanek · 29 Jul 2012

Try again with lin frequency scaling (Bode Plot properties).

infuzion · 29 Jul 2012

martinvicanek wrote: Try again with lin frequency scaling (Bode Plot properties).

Has the same zig-zag to me, just moved over to the right.

martinvicanek · 29 Jul 2012

@infuzion, the zig-zag is because the display is confined to a range between -180 deg and 180 deg. However note that there is an ambiguity as you can always add multiples of 360 deg to the phase without changing anything. If you unwind the graph in this way it will actually become linear. (That holds for the passband. In the stopband things are a bit more complicated.)

Warmonger · 30 Jul 2012

Dude, but the range of zig-zag is neither contstant nor full. It's also full of discontinuities.

infuzion · 30 Jul 2012

I did a bit of re-arranging so you can change the Phase Range on the front panel Preferences for the Bode Plot. I did not see any real difference, aside from the obvious ratio-squishing. Maybe I did something wrong?

Also, when the filter's frequency input is changed, seems the initial phase change does not change.

martinvicanek · 30 Jul 2012

range of zig-zag is neither contstant nor full. It's also full of discontinuities.

There are zeros in the stop band that result in a 180 deg discontinuity each. It's quite common if you look here, here, or here (scroll down for the last two).

Maybe I did something wrong?

No, it's just that the GraphFF module will always output a phase value between -180 and 180.

Also, when the filter's frequency input is changed, seems the initial phase change does not change.

Yep, that's because the slope represents the group delay which, in the example posted, is 16 samples regardless of the cutoff.

Warmonger · 30 Jul 2012

Ok, I see the problem. Bode plot does not draw actual phase of complex signal, just the phase of magnitude.

infuzion · 30 Jul 2012

Warmonger wrote: Ok, I see the problem. Bode plot does not draw actual phase of complex signal, just the phase of magnitude.

OH! Martin, why didn't you say that in the first place! 😉

infuzion · 6 Aug 2012

martinvicanek wrote: infuzion wrote: BTW, I would pre-calculate n - 16; though not much of a speed increase on its own due to the hop(), it has additional benefits, like buffering, helping with ASM compiler optimization, etc.

Oh, sure, there is plenty of room for optimization! For instance it is a complete waste to use sin1() and cos1(), much better to compute recursively.

Do you mean like this? //nX = 16-X; w = (1+ (cos1(nX /32))) / nX; float n1 = 15; float w1 = 0.0012809813064513047; float n2 = 14; float w2 = 0.005437176249193804; float n3 = 13; float w3 = 0.012963875976727282; float n4 = 12; float w4 = 0.024407768234454376; float n5 = 11; float w5 = 0.0404027060891271; float n6 = 10; float w6 = 0.06173165676349103; float n7 = 9; float w7 = 0.08943440866487463; float n8 = 8; float w8 = 0.125; float n9 = 7; float w9 = 0.1707271888594469; float n10 = 6; float w10 = 0.23044723872751496; float n11 = 5; float w11 = 0.31111404660392045; float n12 = 4; float w12 = 0.42677669529663687; float n13 = 3; float w13 = 0.610489870767515; float n14 = 2; float w14 = 0.9619397662556434; float n15 = 1; float w15 = 1.9807852804032304;

hop(4096) // filter coefficients for windowed sinc { //green fh = 0.5*f; c1 = w1 *sin1(n1 *fh); c31 = c1 ; c2 = w2 *sin1(n2 *fh); c30 = c2 ; c3 = w3 *sin1(n3 *fh); c29 = c3 ; c4 = w4 *sin1(n4 *fh); c28 = c4 ; c5 = w5 *sin1(n5 *fh); c27 = c5 ; c6 = w6 *sin1(n6 *fh); c26 = c6 ; c7 = w7 *sin1(n7 *fh); c25 = c7 ; c8 = w8 *sin1(n8 *fh); c24 = c8 ; c9 = w9 *sin1(n9 *fh); c23 = c9 ; c10 = w10 *sin1(n10 *fh); c22 = c10; c11 = w11 *sin1(n11 *fh); c21 = c11; c12 = w12 *sin1(n12 *fh); c20 = c12; c13 = w13 *sin1(n13 *fh); c19 = c13; c14 = w14 *sin1(n14 *fh); c18 = c14; c15 = w15 *sin1(n15 *fh); c17 = c15; c16 = fh*12.56637061435917295385; //n=16;w=2;c16=w*fh*2*3.14159; treat separately to avoid div by 0 norm = c1+c1+c3+c4+c5+c6; norm = norm+c7+c8+c9+c10+c11+c12; norm = norm+c13+c14+c15+0.5*c16; norm = 0.5/norm; } Saved like 10 whole cycles on my i7. Could have gotten that just by increasing the hop(). 😮

trogluddite · 6 Aug 2012

The two biggest savings surely would be... - Co-efficients are symmetrical - so the 'mirror image' coeff's do not need separate variables (memory) - and the maths could then be factorised to lose nearly half of the multiplies. (out = (in1 + in31) * c1; etc...) - Hate to say it, but... ...this one is ripe for ASM conversion - the compiler is ridiculously wasteful when repeating operations on the same variable - every 'out+out+...' does a totally redundant read and write when, instead, the running total could be kept in an xmm register. For simplicity, I'd just copy and paste the hopped part. it could be optimised but at hop(4096) it is probably not worth the considerable effort involved - however, the main iteration body would make a good 'entry level' optimisation exercise as it is all straightforward SSE maths. (see the Wiki for a good tutorial on doing this).

martinvicanek · 6 Aug 2012

Totally agree with you, Trog! Infuzion, your idea to precompute the cosines for the window function is even better than what I had in mind! However, for the remaining sines one could use a recursive scheme (calculate the current sine easily from the previous one, like here). In fact one should also precompute the 1/n and factor them into the window function. Edit: Oh, I see you have done that already. Yet another improvement would be to use a circular buffer instead of fire brigade shifting samples. Indexed variables are inefficient in code, but in ASM Trog has some tricks up his sleeve IIRC.

trogluddite · 7 Aug 2012

martinvicanek wrote: Indexed variables are inefficient in code, but in ASM Trog has some tricks up his sleeve IIRC.

Well, indexed buffers can certainly be more optimal that the stock code - but the bucket brigade is usually still much more efficient in a case like this where you need to read multiple values from the buffer at each sample. No matter how well the indexing is optimised, you always have to split the SSE channels when using indirect addressing, and then read memory in four groups of 32bits via a temp variable - rather than the nice chunky 128bit "all-in-one" SSE reads direct to a register.

infuzion · 7 Aug 2012

trogluddite wrote: For simplicity, I'd just copy and paste the hopped part. it could be optimised but at hop(4096) it is probably not worth the considerable effort involved - however, the main iteration body would make a good 'entry level' optimisation exercise as it is all straightforward SSE maths. (see the Wiki for a good tutorial on doing this).

I did this at first, but decided to optimize inside the hop() because: 1) In case I wind up dropping the hop() because attaching a LFO/env to the cutoff 2) Remove the redundant symmetrical variables (thanks) 2) Save file/RAM space from a smaller ASM coding But it seems slightly slower by 2 cycles or so afterwards for some reason. I might have to do ninja opcode interpolations or something.

After all the ASM optimizations, I got the cycles on my i7 from ~350 down to ~82; one of my better % drops. Could drop it more of course...

martinvicanek wrote: Infuzion, your idea to precompute the cosines for the window function is even better than what I had in mind!.. ...Yet another improvement would be to use a circular buffer instead of fire brigade shifting samples. Indexed variables are inefficient in code, but in ASM Trog has some tricks up his sleeve IIRC.

Thanks! Best to always look for static values for reduction.

Not sure what you mean by "circular buffer"; do you mean an array[x]? -- Warning: The Transfer Function (Frequency Domain) module seems to lock up opening files; I suggest changing the Tick 25 division to 5 so there are only 5 updates/second or so.

martinvicanek · 7 Aug 2012

infuzion wrote: by "circular buffer"; do you mean an array[x]?

Yes but see Trogs comment.

infuzion wrote: Warning: The Transfer Function (Frequency Domain) module seems to lock up opening files; I suggest changing the Tick 25 division to 5 so there are only 5 updates/second or so.

Thanks!

trogluddite · 7 Aug 2012

infuzion wrote: I got the cycles on my i7 from ~350 down to ~82

Result! 🙂

infuzion wrote: 1) In case I wind up dropping the hop() because attaching a LFO/env to the cutoff

Ah yes, quite right - hop(4096) could be zippy as a punk rocker's leathers!

infuzion wrote: But it seems slightly slower by 2 cycles or so afterwards for some reason. I might have to do ninja opcode interpolations

Hard to say without seeing your code, but before donning your black balaclava, consider that this whole routine is chock-a-block with memory reads/writes - so fine tuning the cache usage may be more critical than optimising for out-of-order processing (which will still help to free up the instruction pipeline during what cache misses remain). One suggestion that's quick to try out... - Due to the factorising, and the chain of "In{N} = In{N-1}", the "In" variables often occur in the code in reverse order - so it may be worth declaring the floats in reverse order in the code header, so that the memory addresses are read more sequentially. No guarantees, but I have found some good savings this way in some buffer algorithms, due to fewer cache misses (even one miss can mean a stall for dozens of cycles).

PS) Some great (and fairly easy reading) cache optimising tips HERE. I bookmarked this when I was re-building my SoopaLoopa (dozens of 30sec buffers), and following the advice (and links) in there got me a 30% CPU reduction before I even started slaughtering opcodes.

infuzion · 8 Aug 2012

trogluddite wrote: ...fine tuning the cache usage may be more critical than optimising for out-of-order processing (which will still help to free up the instruction pipeline during what cache misses remain)... - Due to the factorising, and the chain of "In{N} = In{N-1}", the "In" variables often occur in the code in reverse order - so it may be worth declaring the floats in reverse order in the code header, so that the memory addresses are read more sequentially...

Thanks for the tip! I knew recent-accessed is better than long-ago-accessed, but I didn't think about memory order too much. All I really can do is reverse the order for the inX variable for the actual filter, and move them to the end of declarations (unless they should be listed first due to the hop()?). How can I test the speed difference please? Signal Analyzer seems to be of little help...

trogluddite · 8 Aug 2012

infuzion wrote: (unless they should be listed first due to the hop()?

I think you're probably right to declare the variables used only in the hop first - just because it is earlier in the code. Otherwise, I don't yet have enough experience with this myself to offer any hard-and-fast predictions. I would guess that there's little reason to treat hopped code much differently to regular code - at SR=44.1kHz vs. CPU=??GHz, even normal code leaves plenty of time for the variables to be completely flushed from the cache between samples - e.g. by other plugins, the VST host, OS etc.. So I would think that ordering things for best performance of a single unhopped iteration is the most important factor - but with all sorts of modern CPU pre-fetch trickery going on, I could well be wrong.

infuzion wrote: How can I test the speed difference please? Signal Analyzer seems to be of little help...

A tricky one - even if the analyser were more accurate, "cache effects" are highly dependent on what else is going on with the memory (other processes, or even other modules within the same schematic). I have a few times made modules that appeared to be less efficient when tested in isolation, only to find that they were much more effective when used in multiples, or as a plugin within a host (and other time the exact reverse!) So, although it's more tedious, I tend to use "real world" tests - such as running large multiples of the module to get a nice "average" CPU meter reading. and ultimately, using Reaper's "per track" CPU metering to test a chain of VST exports. For example, the SoopaLoopa that I used as an example was quite simple to test, as there were so many identical buffers within the design - running 10 instances in Reaper gave me nice steady CPU loads well into double figures, making for very simple comparisons. And, of course, just like out-of-order optimisation, so much will depend on the precise architecture of the PC being used for testing - I doubt that there is any kind of test that will translate well between different CPUs.

Dell · 8 Aug 2012

If you want to test the speed, the best thing to use is a tone burst generator. It will allow you to set the decay level based on your requirements.

This should be able to get the job done. http://trueaudio.com/rta_abt1.htm

infuzion · 8 Aug 2012

trogluddite wrote: infuzion wrote: How can I test the speed difference please? Signal Analyzer seems to be of little help...

A tricky one ...I tend to use "real world" tests - such as running large multiples of the module to get a nice "average" CPU meter reading. and ultimately, using Reaper's "per track" CPU metering to test a chain of VST exports.

I think that is the best bet. cheers