Biquad coefficients using Magnitude (or Phase) Invariance Mapping Method

[juha_p]

I'm trying MIM (Magnitude Invariance Method) and PIM (Phase Inveriance Method) for to improve biquad LPF response at low sampling rates. I'm not sure if these methods can be used for this purpose so I'm looking some help and examples of usage if available. Paper (Google finds the pdf where both methods are described (soar.wichita.edu -link).

[keithwood]

I always found that any algorithm that relies on cepstral methods does not work if the analog transfer function is zero at any point (such as the center frequency of a notch or zero frequency for a high pass) because it then has to take the log of zero (or a very small number). If anyone has a way round that I'd be interested to know.

[juha_p]

Instead of controls (as in paper) my aim is to try these methods for to improve biquad LPF 'output' at low sampling rates. Since my knowledge level regarding DSP and involved math is quite low (i.e. that paper tells me not much), finding some source code (C/C++) would help me to start do some experiments.

[keithwood]

If it's a lowpass that has a more analog magnitude response, you want to search for the Massberg lowpass. This approach models a 1st or 2nd order low pass via an algorithm that choose appropriate frequency and Q and gain for a high shelf. It's probably the best you're going to get if you're short on time/knowledge - it's a lot better than using the bilinear transform.

The paper is available from AES under the title "Digital Low-Pass Filter Design with Analog-Matched Magnitude Response" which from memory I think has some matlab code with it. Also Will Pirkle's book Designing Audio Effect Plug-ins in C++ with Digital Audio Signal Processing has the design equations in it which are simply laid out.

[juha_p]

Thanks for suggestions. I have the Massberg already implemented but would like to try other available methods as well. What comes to those DSP books, my broblem with them is the foreign language (even papers/books written by finns are written mostly in english).

I have not yet tried the Matlab example found in MIM/PIM paper. To run the example code in Octave for LPF what 'general' changes are needed to do to get valid results (Matlab vs Octave synthax and input parameters)?

[juha_p]

OK, I tried the Matlab listings found in paper in Octave for MIM method. No errors reported but the plots are OK only for the analog model.

Never mind ... found the reason for issue and now plots are OK.

[keithwood]

Glad you're getting somewhere. I use both matlab and octave and I've found that scripts written in one will run in the other unless you're using a built in function that only exists in one of them. I can't see anything in the MIM paper that looks like it will cause octave any problems.

[juha_p]

Glad you're getting somewhere. I use both matlab and octave and I've found that scripts written in one will run in the other unless you're using a built in function that only exists in one of them. I can't see anything in the MIM paper that looks like it will cause octave any problems.

Yes, there were zero errors reported by octave. Required sampling time ('Ts') was set to 1 in paper and it gave different results from the ones put in paper:
http://i66.tinypic.com/33jpx7b.png
By using lower values (as like 1/10 ... 1/50) gives 'better' results:
http://i67.tinypic.com/nsug0.png

Latter plots looks similar with the figures (4.12-4.15) in paper but not equal so something still needs to be looked at (the paper is copy protected so re-writing the code might lead to typos).

I haven't got the PIM part of paper working yet. There's variable 'i' in code line

Code: Select all

PhaseJ = Phase * i

which of I have no glue where it comes from. Tried few values for 'i' but nothing changed so there might be something wrong with the parameters (script is originally coded for MIM).

In original Matlab script there are two sets of coefficients given ... for 'plant' and 'analog' filter. How these should be set for LPF filter ... what 'plant' means there ?
Also, what should be done to get the frequency units changed from rad/s to Hz units?

[keithwood]

The variable "i" is a built in variable in matlab and octave and is the square root of minus one (imaginary one).

In a paper like this a "plant" is a filter of known coefficients (typically in digital domain) which different models can be tested against to see how closely they can approximate it, then relative errors can be determined.

[juha_p]

The variable "i" is a built in variable in matlab and octave and is the square root of minus one (imaginary one).

In a paper like this a "plant" is a filter of known coefficients (typically in digital domain) which different models can be tested against to see how closely they can approximate it, then relative errors can be determined.

Thanks for clearing these.

OK, some progress.
Tried the MIM method for 2nd order
- LP filter (fs=44100/f0=15000/Q=0.707) and
- peak filter (fs=44100/f0=10000/Q=0.707/db=6).
Here are the plots:
http://oi68.tinypic.com/34qlkzb.jpg
http://oi67.tinypic.com/2ymf2me.jpg

Magnitude responses looks good otherwise but the magnitude and/or frequency areas are not what I would like to have (any problem ?). How can I 'move' the filters to work at normal range (dB/Hz)? I'm using Octave here. Also the 2nd order peak filter (MIM) has little more gain than set +6db ... what could cause this?

[juha_p]

There are few issues which needs to be solved:

• 'move' the filters to work at normal range (0-fs/2)? (SOLVED)
• fc frequency effects onto opposite direction -> now wrong since when fc is set to 10kHz resulting filter coefficients plot it at ~5kHz, 15kHz -> ~7.5kHz (SOLVED)
• gain normalization (Lp filter) (SOLVED)
• gain matching (2nd order peak filter is off a bit as seen in plot)

Any thoughts how to solve these issues?

Example code (LPF):

Code: Select all

% Analog model for LP filter:
fs = 44100;
fc = 10000;
Q = 0.707;

w0 = 2 * pi * fc;
b2 = w0^2;
a2 = [1 w0/Q w0^2];
Ds = tf(a2, b2);

% Calculate coefficients using MIM
[Dz] = c2dn(Ds, fs, 'mim', 2, 4096*256);

bodemag(Ds, 'r', Dz, 'b');

[juha_p]

LPF result looks good as seen in plot (frequency range upto 23873Hz):


I have not plotted the error yet.
Coefficients for 2nd order MIM LP filter showen in plot:

Code: Select all

# Created by Octave 4.0.3, Sun Jul 24 12:34:52 2016 GMT <unknown@Juha-PC>
# name: a2
# type: cell
# rows: 1
# columns: 1
# name: <cell-element>
# type: matrix
# rows: 1
# columns: 3
 0.437417389844334 0.1220119196463959 -0.04071587087954469

# name: b2
# type: cell
# rows: 1
# columns: 1
# name: <cell-element>
# type: matrix
# rows: 1
# columns: 3
 1 -0.6951293456643256 0.2138427842755107

# name: scalarDs
# type: scalar
1

# name: scalarDz2
# type: scalar
1

So, it looks like MIM can be used for to improve 2nd order LPF magnitude response.

[juha_p]

I always found that any algorithm that relies on cepstral methods does not work if the analog transfer function is zero at any point (such as the center frequency of a notch or zero frequency for a high pass) because it then has to take the log of zero (or a very small number). If anyone has a way round that I'd be interested to know.

I'm noob with math and math programming so I have to ask if algorithm for calculating log (x) could be used for to round that issue?
In case of MIM/PIM Matlab coding there is these code lines

Code: Select all

R = [R1 R2];  %% magnitude squared response
R = R';
h = log(R); % R^(k)

,... how would that be handled if that polynomial approximation is used instead of log()?

[keithwood]

I meant that log(0) = -infinity so you have to take care, for low pass should be fine though.

[juha_p]

I meant that log(0) = -infinity so you have to take care, for low pass should be fine though.

Octave crashes now by HP and BP (probably because of this log() issue or something related to it and dcgain() function gives warnings as well). Looks like 'cookbook' filter types can be calculated using MIM except those two. AP gets calculated but phase response is wrong and magnitude response is wrong for low order Notch filter when value of fc is low.
I was thinking to try if that polynomial approximation way could be useful in getting build those 2 filters but would need some help coding it for Octave use.