What could be the reason for this 'error'

[juha_p]

Hi!

Here's MIM (Magnitude Invariance Method) paper (press the "View PDF" to get the paper/matlab code in hands) and here's a plot showing the issue I'm trying to solve:



Magnitude responses for one pole LP filter at few frequency locations.

As plot shows, there's some issue in MIM's function c2dn() implementation, which 'moves' the fc point especially at low fc values. IIRC, same type of issue happens when changing the fs.

Analog prototype LPF is calculated as:

Code: Select all

w0 = 2*pi*fc;
Analogb = 1;
Analoga = [1 w0];
Analog = tf(Analogb, Analoga);
Analog = Analog/dcgain(Analog);

which is then given as an input for to calculate the MIM filter:

Code: Select all

fs=44100;
mimpim = 'mim';
samplingtime = 1/fs;
numofsamples = 2048;
[Dz1] = c2dn(Analog, samplingtime, mimpim, 1, numofsamples);
%[a2, b2, T] = tfdata(Dz1);   %coeffs

where c2dn() is the function from MIM paper.

Plotting:

Code: Select all

fs2 = fs/2;
nf = logspace(0, 5, fs2);
[mag, pha] = bode(Analog,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'g', 'linewidth', 1.5);
axis([10 fs2 -50 1]);
hold on;[mag, pha] = bode(Dz1,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'k', 'linewidth', 1.5, 'linestyle', '--');grid on;
title('Various TF (LPF 1)');
legend('Analog', 'MIM','location', 'southwest');
xlabel('Hz');ylabel('dB');

So far I found out that the fc error is quite linear between certain points (I picked -3dB points for few fc and then calculated the trend f(x) (R² was 0.99989...)) so, one could probably see the source for "error" by looking the matlab codes for MIM function c2dn().

Here's plot showing the effect of my correction f(x):



Correction formula

Code: Select all

f(x) = 0.9704384746x - 51.6991952722

works well for fc in range 60Hz-Nyqvist (fs=44.1kHz). (Note: correction showen in plot for fc=10Hz is done using another formula so you can omit it as a result).

Any suggestions regarding the source for this "error"?

[mystran]

I looked at the paper and it's too long for me to read properly right now.. but I just wanted to note that..

If you have a frequency normalised analog transfer function like 1/(s+1) then you can tune this to any frequency w0 by replacing s by s/w0 (and s^2 by (s/w0)^2 and so on). The nice thing about this substitution is that it doesn't change the DC magnitude (or well, any magnitude, it just scales the frequency). It is easy to verify that 1/(s/w0+1) = w0/(s+w0), but it works for arbitrary filters of arbitrary orders, with or without fancy zeroes.

[juha_p]

Thanks for the tip.

Decreasing the value of numofsamples parameter improves the magnitude response at low fc (and vice versa) ? Could the issue be related to those fregs() calls in mim block or to cepstral processing block?

Note: When you build HPF (low fc and fs) using MIM, its needed to set the numofsaples parameter value to its upper limits to get good magnitude response (even when filter order is > 2).

EDIT: I quess the functionality behind the issue can be found from cepstral processing block ... there's a line with remark: "% minimum phase sequence r^mn". When this line is in comments magnitude plots looks good (if the numofsamples parameter has suitable (big) value):



Problem is that commenting the line effects directly to the aliasing issue this type of filter has

Any thoughts if this (r^mn thingy) can be improved somehow?

[juha_p]

When the resulting LPF order is 3 or more this issue in discussion is not present anymore:


(numofsamples=4096, fs=44.1kHz, fc=1...N Hz)

--> could it just be so that the fitting method used there in function c2dn() isn't suitable for low order filter?

[DaveClark]

Assuming you have not already fixed this problem:

If I understand your plots correclty, they show markedly increasing problems with lower cutoff frequencies, yet this is where all the author's results come from --- in fact, far lower. If you have not already done so, I recommend that you first try to duplicate the author's results with his/her data to ensure you have not introduced some typo or are misunderstanding some units somewhere. You've probably got something way the heck off that easily reveals itself in the low-frequency limit.

Good luck.

[juha_p]

Assuming you have not already fixed this problem:

If I understand your plots correclty, they show markedly increasing problems with lower cutoff frequencies, yet this is where all the author's results come from --- in fact, far lower. If you have not already done so, I recommend that you first try to duplicate the author's results with his/her data to ensure you have not introduced some typo or are misunderstanding some units somewhere. You've probably got something way the heck off that easily reveals itself in the low-frequency limit.

Good luck.

Yes, I have not fixed it ... just rounded it.
IMO, you be able to see the issue in thesis plots as well. Lower fc ... yes, but the fs is lower there as well (I think it's question of fc:fs:samplebuffer).

Actually I did run few of the examples given in thesis (also got some help here for to solve couple of them (IIRC, wrong data was given in thesis text)).

I have checked few times already that my Octave source is 1:1 with the Matlab source given in thesis. I can't say much regarding the Matlab implementation author did (dunno all what's happening there) nor if there are differences between Matlab and Octave in functions used there but, as I told already, code works perfectly when filter order is set to 3 or higher and databuffer is big enough for selected fc:fs.

[AnalogGuy1]

Ask the author? https://www.linkedin.com/in/prathamesh- ... r-5200314/

Someone suggested a while ago when I was having trouble creating a patch for a Car's song to just ask Greg Hawkes. Friended him on Facebook and asked; within a few hours I had the answer. Since then I've done it a number of times when I've had trouble understanding something in a book.

[juha_p]

Ask the author?
...

I quess chapter 5 there in thesis kind of explains some issues and needed future work.
- higher order filter is needed when sampling time is decreased (he uses 1 second in examples ... I have it usually set to 1/44100 seconds to 1/384000 seconds
- something related to periodic property H(e^jw) - ?
- cusp in magnitude observed around ±π

I have not found future work papers for MIM/PIM methods.