How to create peak/shelf filters, that are symmetric close to nyquist?

Ok I made it work, but the shape isn't still perfect, now I veryfied in
Octave.

Basically none of the Q definitions work well here. In 44k sampling rate,
it mostly looks ok, when the selected frequency is below say 10k. Above it
the shape starts getting "thin", maybe not as normal filter, but still a
lot. Also when Q < 0.4, it starts doing really weird stuff. For example
when using Q = 0.05 and moving the frequency accross the spectrum, the
"thickness" of the peak filter is getting up and down, more or less waving.
It looks to me like the approximation of the requested gain in Nyquist may
not be right. Any ideas?

jungledmnc

On 11/10/11 8:28 PM, jungledmnc wrote:
> Ok I made it work, but the shape isn't still perfect, now I veryfied in
> Octave.
>

lessee... what did i say?

On 11/3/11 11:05 PM, robert bristow-johnson wrote:
 >
 > but there is no way to make it perfect because the slope of the gain
 > function at Nyquist must necessarily be 0 for a digital filter (and it
 > isn't for the analog prototype).
 >



> Basically none of the Q definitions work well here.

how are you defining Q?  in terms of apparent BW?  like this:

     1/Q = 2*sinh( ln(2)/2 * BW * w0/sin(w0) )   BW in octaves

or by some other definition?

is scaling up bandwidth (however you're defining it) by w0/sin(w0) 
sufficient?  maybe not so when w0 gets really, really close to Nyquist?

> In 44k sampling rate,
> it mostly looks ok, when the selected frequency is below say 10k. Above it
> the shape starts getting "thin", maybe not as normal filter, but still a
> lot. Also when Q < 0.4, it starts doing really weird stuff.

consider what the corresponding BW is.  what happens if the bandedge on 
the right exceeds Nyquist?  would you expect "normal"?

> For example when using Q = 0.05

how wide is BW when Q = 0.05?

> and moving the frequency across the spectrum, the
> "thickness" of the peak filter is getting up and down, more or less waving.
> It looks to me like the approximation of the requested gain in Nyquist may
> not be right. Any ideas?

maybe set it to be a little higher gain (assuming boost) than what 
Orfanidis does.

jungle, you're pushing the edges a little here.  expect wierdisms or 
change Fs to 96 kHz (from which the regular old cookbook works fine even 
at around 20 kHz).

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

On 11/10/11 8:23 AM, Rune Allnor wrote:
>
> Matlab&  clones usually take a lot of beating in this
> forum, so for once I will praise it:
>
> This is exactly the kind of thing matlab (or octave or
> scilab or...) was invented for!

fine, so then what's The MathWorks excuse for doing such a bad job of it?

> You, the user, reads some document about some ethod
> and want to try it. Instead of diving deep into C or
> Fortran or assembler or whatever straight away, you
> first use matlab (or octave or scilab or...) to check
> that you get the computations right.
>

before i had access to MATLAB (i have never purchased it), i wrote a 
quickie little C program for my Mac that would plot the graph of 
whatever real function i could define in C.  i would have to recompile 
the "f.c" file, link, and run.  this worked very nicely with Lightspeed 
C (later to become THINK C, still later Semantec C, then death at the 
hands of Code Warrior), just open the project, type in an equation and 
run.  i had FFT (it was before the FFTW days) and some audio in and out 
utilities (AIFF and WAV).  the compile time was negligible (compared to 
0 for an interpreted language like MATLAB).  execution time was very 
fast.  and the DC bin is X[0].

i do not denigrate the need or purpose of MATLAB.  i am fully 100% 
critical of their indexing scheme; mostly that the index base is 
hard-wired (where even the array dimensions are not hard-wired).  if 
MATLAB has a function called "reshape()", why can't it have a function 
called "reorigin()" or "rebase()"?  there is no good reason why not and 
the fact that after decades they *still* put the DC bin into X(1) means 
they haven't gotten it.  that cannot ever be considered "right".

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

On Nov 10, 6:41 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> ...
>
> i do not denigrate the need or purpose of MATLAB.  i am fully 100%
> critical of their indexing scheme; mostly that the index base is
> hard-wired (where even the array dimensions are not hard-wired).

Actually, the dimensions are hardwired. The first dimension is 1.
You've re-whined this so many times your abstractions are blurring.

>  if
> MATLAB has a function called "reshape()", why can't it have a function
> called "reorigin()" or "rebase()"?

It could, but then there would be users writing C in Matlab as well as
Fortan in Matlab, when the point of a language like Matlab is to avoid
low level code whenever posible, not facilitate it.

> there is no good reason why not and
> the fact that after decades they *still* put the DC bin into X(1) means
> they haven't gotten it.  that cannot ever be considered "right".
>

It has been and always will be right for Matlab's primary marketing
audience, unless you can convince all authors from Wikipedia to Golub
and van Loan to teach 'Matrix Computations' in zero base index. The
Mathworks also has the right idea that they wouldn't be getting any
more subscription payments from whiney old C programmers on comp.dsp
if they had done things differently.

Dale B. Dalrymple

On 11/11/11 2:30 AM, dbd wrote:
> On Nov 10, 6:41 pm, robert bristow-johnson<r...@audioimagination.com>
> wrote:
>> ...
>>
>> i do not denigrate the need or purpose of MATLAB.  i am fully 100%
>> critical of their indexing scheme; mostly that the index base is
>> hard-wired (where even the array dimensions are not hard-wired).
>
> Actually, the dimensions are hardwired.

actually, no they're not.  that's what the reshape() function changes. 
your array can be changed from 2x6 to 3x4 or to 12x1.  these limits are 
checked by MATLAB every single time that arrays (or "matrices" if you 
want) are added or multiplied or when individual elements are accessed.

so the upper limit can be assigned, why not the lower limit (besides the 
fact that MATLAB deigns not to provide this functionality).

> The first dimension is 1.
> You've re-whined this so many times your abstractions are blurring.
>

it should have been obvious, when they were developing the Signal 
Processing Toolbox (or the fft() function) that when DC goes into X(1), 
when the results of the fft() or of conv() come out *wrong* (at least 
differently than any textbook and the published lit), *something* should 
have been nagging the consciousness of Cleve Moler that something wasn't 
quite right.

and they could have fixed this and be perfectly backward compatible.

Dale, if this argument was at comp.soft-sys.matlab, you might have 
friends, but the comp.dspers have long ago figured this out.

>>   if
>> MATLAB has a function called "reshape()", why can't it have a function
>> called "reorigin()" or "rebase()"?
>
> It could, but then there would be users writing C in Matlab as well as
> Fortan in Matlab, when the point of a language like Matlab is to avoid
> low level code whenever posible, not facilitate it.

oh, so when using MATLAB, it "facilitates" us having to subtract 1 from 
the index when we start doing math on it, and when our math is done and 
we come up with another index, we have to add 1 before using it.  such 
"facility", such "utility".

>
>> there is no good reason why not and
>> the fact that after decades they *still* put the DC bin into X(1) means
>> they haven't gotten it.  that cannot ever be considered "right".
>>
>
> It has been and always will be right for Matlab's primary marketing
> audience,

really?  i wonder if more EEs doing signal processing use MATLAB than 
any other single group.  the Signal Processing Toolbox was the very 
first toolbox that they came out with (in the early 1990s).  we may not 
be the majority of users, but i'll bet that EEs are the plurality.

> unless you can convince all authors from Wikipedia to Golub
> and van Loan to teach 'Matrix Computations' in zero base index.

you still haven't gotten it, Dale.  can you spell "backward compatible"?

> The
> Mathworks also has the right idea that they wouldn't be getting any
> more subscription payments from whiney old C programmers on comp.dsp
> if they had done things differently.

it's a losing argument if merit is the basis.  but if corporate inertia 
is the basis, you win.


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

>> Basically none of the Q definitions work well here.
>
>how are you defining Q?  in terms of apparent BW?  like this:
>
>     1/Q = 2*sinh( ln(2)/2 * BW * w0/sin(w0) )   BW in octaves
>
>or by some other definition?
>
>is scaling up bandwidth (however you're defining it) by w0/sin(w0) 
>sufficient?  maybe not so when w0 gets really, really close to Nyquist?

I think the definition of Q is actually the problem. To be honest, I know
nothing about this, so I'm just messing with definitions other people
invented. Btw. is there some paper about this? Just something to understand
the basics. 
Anyway I tried these:

const double Dw = Q * w0 / msin(w0);
const double Dw = 2 * w0 * msinh( (msin(w0)/w0) * masinh (1/(2*Q)));
const double Dw = masinh(0.5/Q) * (2/mln(2.0)) * w0 / msin(w0);
const double Dw = msin(w0) / (2*Q);

None of them work well. Changing frequency changes the bandwidth too,
except for the last one, where it works more or less fine, but when the
frequency gets close to nyquist, the shape gets "thin" again.
But in most cases I can change the Q to make it look good again, so I think
the whole trouble here is to find the good Dw definition.
Any ideas?


>consider what the corresponding BW is.  what happens if the bandedge on 
>the right exceeds Nyquist?  would you expect "normal"?

Ideally the bandwidth would still be the same, so that in log view the
shape doesn't change much when changing frequency. It is possible, several
implementations use that.


>> For example when using Q = 0.05
>
>how wide is BW when Q = 0.05?

With Q = 0.05, SR 44k, the filter affect the whole relevant range from 20
to 20k with any frequency. It's kinda extreme, but with does the same thing
with lower Q as well, it's just less visible.


>jungle, you're pushing the edges a little here.  expect wierdisms or 
>change Fs to 96 kHz (from which the regular old cookbook works fine even 
>at around 20 kHz).

I know, but the whole idea is to make it good even with lower sampling
rates such as 44k. And I know it is possible, just don't know how :(.

jungledmnc

On 11/11/11 3:43 PM, jungledmnc wrote:
...
>> jungle, you're pushing the edges a little here.  expect wierdisms or
>> change Fs to 96 kHz (from which the regular old cookbook works fine even
>> at around 20 kHz).
>
> I know, but the whole idea is to make it good even with lower sampling
> rates such as 44k. And I know it is possible, just don't know how :(.
>

*what* is it that you know is possible, and how or why do you know that 
it's possible?


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

>
>*what* is it that you know is possible, and how or why do you know that 
>it's possible?
>
>

That there are several software plugins that implement it without
oversampling (for instance Ozone 5, IIEQ...). It's not exactly perfect, but
it almost perfectly keeps the peak filter shape across the spectrum.

jungledmnc

On 11/11/11 9:26 PM, jungledmnc wrote:
>>
>> *what* is it that you know is possible, and how or why do you know that
>> it's possible?
>>
>>
>
> That there are several software plugins that implement it without
> oversampling (for instance Ozone 5, IIEQ...). It's not exactly perfect, but
> it almost perfectly keeps the peak filter shape across the spectrum.
>

well, maybe they're using a higher order IIR, or maybe they're using an 
FIR (and the FIR need not be linear phase), or maybe they're using a 
combination like a biquad with an FIR "fixing it", but if it's a basic 
2nd-order IIR (a biquad either in straight-forward Direct 1 or in the 
state-variable or some other topology), there are only 5 numbers you can 
play around with.

and no matter how it's done (if it ain't oversampled), the *slope* of 
the gain function must be 0 at Nyquist.  the curve has to start leveling 
off at frequencies just below Nyquist.  you cannot beat that without 
moving Nyquist to a higher frequency, which is oversampling.

so, i'm telling you, jungle, there are limits with what you can do with 
5 knobs and the standard biquad.  this (and the Orfanidis thingie) came 
up several times in different contexts at the last AES convention.  and 
i was sorta surprized.  one reason (for being surprized) is, for the 
5-coefficient biquad, i thought that Knud Christianson of TC electronic 
has done in the previous decade a generalization of the 2nd-order 
peaking and shelving filters (with there are 4 degrees of freedom not 
counting constant gain).  that's a total of 5 degrees of freedom; 5 
coefficients, sounds like a totally solved mathematical problem to me.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

>well, maybe they're using a higher order IIR, or maybe they're using an 
>FIR (and the FIR need not be linear phase), or maybe they're using a 
>combination like a biquad with an FIR "fixing it", but if it's a basic 
>2nd-order IIR (a biquad either in straight-forward Direct 1 or in the 
>state-variable or some other topology), there are only 5 numbers you can 
>play around with.
>
>and no matter how it's done (if it ain't oversampled), the *slope* of 
>the gain function must be 0 at Nyquist.  the curve has to start leveling 
>off at frequencies just below Nyquist.  you cannot beat that without 
>moving Nyquist to a higher frequency, which is oversampling.
>

Yes, actually the filters in all the mentioned software all have "0 slope"
at Nyquist and I'm ok with that. And the Orfanidis filter looks
suspiciously similar! The only trouble is how the filter reacts to
different frequency/Q combinations. With all frequencies (at least so far)
I can find a Q (or rather Dw), which makes it look correct. So I think all
it needs is a good Dw definition. Unfortunately that's exactly what I don't
know anything about :(...

jungledmnc