Biquad coefficients using Magnitude (or Phase) Invariance Mapping Method
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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).
Last edited by juha_p on Sat Jul 23, 2016 7:17 am, edited 2 times in total.
- KVRist
- 91 posts since 24 Dec, 2015 from Bristol, UK
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.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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.
Last edited by juha_p on Tue Aug 16, 2016 7:00 pm, edited 3 times in total.
- KVRist
- 91 posts since 24 Dec, 2015 from Bristol, UK
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.
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.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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)?
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)?
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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.
Never mind ... found the reason for issue and now plots are OK.
Last edited by juha_p on Fri Jul 22, 2016 6:29 am, edited 1 time in total.
- KVRist
- 91 posts since 24 Dec, 2015 from Bristol, UK
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.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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:keithwood wrote: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.
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?
Last edited by juha_p on Fri Jul 22, 2016 6:30 am, edited 2 times in total.
- KVRist
- 91 posts since 24 Dec, 2015 from Bristol, UK
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.
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.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
Thanks for clearing these.keithwood wrote: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.
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?
Last edited by juha_p on Sun Jul 24, 2016 5:44 am, edited 1 time in total.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
There are few issues which needs to be solved:
Example code (LPF):
- '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)
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');
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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:
So, it looks like MIM can be used for to improve 2nd order LPF magnitude response.
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
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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?keithwood wrote: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.
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()?
- KVRist
- 91 posts since 24 Dec, 2015 from Bristol, UK
I meant that log(0) = -infinity so you have to take care, for low pass should be fine though.
-
- KVRian
- Topic Starter
- 834 posts since 21 Feb, 2006 from FI
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.keithwood wrote:I meant that log(0) = -infinity so you have to take care, for low pass should be fine though.
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.