Algebraic minimum-phase variant of sinc function

[Aleksey Vaneev]

camsr, that will be basically the same thing as band-limiting I've described.

In either case, my initial idea was to develop method to be usable for any biquad filter, without using any higher order transformations to avoid 'warping' near Nyquist. Such 'band-limiting' is a straight-forward method, easily understandable.

[camsr]

I wanted to ask what's a good book to read for beginning DSP, like filters. TBH I have no idea what is being said in those graphs with 0 - pi.

[Aleksey Vaneev]

I wanted to ask what's a good book to read for beginning DSP, like filters. TBH I have no idea what is being said in those graphs with 0 - pi.

It's a tough topic for sure... Websearch for 'digital filter design'. However, even a ton of tutorials may not help to understand this topic.

[Aleksey Vaneev]

if you are concerned with alternative stability-checks (other than factorization and checking for poles outside the unit circle), maybe the levinson-durbin recursion is a hot candidate. it converts direct-form coeffs into reflection-coeffs which are esured to produce a stable filter when their magnitude is smaller than 1. ...well, we are drifting far far off topic now...

This is probably the same as decomposing FIR filter of order N into N 1st order filters. Then you can check each 1st order filter for stability easily. I now wonder how such decomposition can be accomplished. (this is probably the same as extracting poles).

[Aleksey Vaneev]

After checking a couple of additional filters I'm leaning to think that stability may be related to the quick phase change at Nyquist. If you consider filter [1 1] - phase jumps quickly at Nyquist (and a bit above) from "pi" to "-pi". Another unstable filter [1 1.5 0.5] exhibits the same behavior. So, if you consider phase response function to be specified in the range 0..2*pi, discountinuity is obviously available.

[Christian Schüler]

A singular value decomposition of the system feedback matrix should do the trick. If all eigenvalues have magnitude < 1, the system should be stable. You should be able to get a feedback matrix trivially if the filter is expressed in state variables.

[Music Engineer]

braindoc, could you plot group delay as well?

here it is:



generated in MatLab's filter visualization tool with the command:
fvtool([1], [1 2])

btw, and the impulse response looks like this:



be sure to recognize the scaling of the y-axis! - the 18th sample is above 25000. that's instability at its best (or worst)

[Music Engineer]

I have no idea what is being said in those graphs with 0 - pi.

it's the frequency response (magnitude and phase) on a normalized frequency axis (normalized with respect to the sample-rate). 0 is dc, pi is Nyquist.

[Aleksey Vaneev]

A singular value decomposition of the system feedback matrix should do the trick. If all eigenvalues have magnitude < 1, the system should be stable. You should be able to get a feedback matrix trivially if the filter is expressed in state variables.

I was actually talking about decomposition of a FIR filter in the first place. So, no feedbacks in it.

[Christian Schüler]

I took the time to read your long post. So its about resampling filter coefficients. That would make for a long number of coefficients then, because the sinc function is infinite.

The IIR coefficients are simply the FIR of the left hand side, you can say, the input filtered by the A coefficients equals the output filtered by the B coefficients. So yeah, both can be resampled the same way.

You will ever only get an approximation to a given filter transfer function at a different sample rate, because at DC and Nyquist, where the conjugate pairs meet, the transfer function must be horizontal. You cannot easily have a clean scissor cut of a filter transfer function in a different sample rate.

[Aleksey Vaneev]

Christian Schüler, where conjugate pairs meet is only a single point (BTW, you do not necessarily have to talk about it in pole/zero design terms - just see how 1st order FIR filter behaves - it always enters Nyquist freq horizontally). But this does not put too strict constraint on how the final response should look (when you divide response of B coeffs by response of A coeffs), there are limits, but as you've said all this does is puts limit to the approximation a bit. - on the other hand A and B coeffs should probably be handled together for best fit.

So, now I have to find method to get rid of a lot of additional coefficients produced after sinc filtering. Probably this step can even be "factorized" leading to a simple algorithm (design filter, sinc-filter, resample, get rid of unnecessary 'sinc-induced' coefficients).

[Aleksey Vaneev]

one more clarification (regarding IIR part): IIR part resembles convolution of a signal with an impulse response which is a deconvolution of the IIR part's coefficients like they were FIR filter taps.

sorry for thinking out loud.

[Aleksey Vaneev]

Getting back to 'resampling' filter design method, it is analytically obvious that such resampling does not cause ANY change in the required final response. When you divide FIR part by IIR part, all freq/phase changes induced by the min-phase sinc cutoff filter (applied to both parts) are cancelled out, and all you have in the end is the response of the original filter, shifted (multiplied) by the resampling coefficient.

Also such cutoff filtering is applicable to IIR part because min-phase filter has all of its zeros inside the unit circle, and so it is applicable in IIR topology (min-phase FIR filters can be used in FIR and IIR parts interchangably).

So, probably this knowledge may help to find a way to 'factorize' or 'algorithmize' such 'resampling' filter design method.

[Aleksey Vaneev]

...The IIR coefficients are simply the FIR of the left hand side...

Took my attention. I think it is some kind of a wrong knowledge.

I've seen somewhere (probably a known 'dspstuff.txt' file - but I'm not sure) something like this: [a0 a1 a2] = [b0 b1 b2].

By definition H(z)=[b0 b1 b2]/[a0 a1 a2] (each component multiplied by z shift, of course).

from this equation you can't say [a0 a1 a2] = [b0 b1 b2], but you should say:
[a0 a1 a2] * H(z) = [b0 b1 b2]

You can't throw H(z) from the equation.

[Aleksey Vaneev]

Hah.. This whole 'coeff' problem can be solved this way:

Take FIR (or IIR part) of the filter, resample it using min-phase sinc function (as I've described). Do the same on dirac spike (FIR filter with a single tap equal to 1.0). Now use this sinc-filtered dirac spike as an IIR part filter. Process that resampled FIR/IIR part of the filter with such filtered dirac spike. This will produce resampled coefficients we are looking for without 'sinc' oscillations (they will be cancelled out). Hopefully resulting filter length will be minimal - but on the other hand, it should become optimal while retaining the original EQ shape.

Viola. (complete same steps with an opposite part of the filter)

note: since we expect the resampled filter to be short, IIR filtering with dirac spike should consist of several iterations only, with a lot of feedback values equal to zero. Beside that we do not have to process the whole 4000 samples of the resulting 'sinc' filtered filter. At least I think so...