Here's MIM (Magnitude Invariance Method) thesis (press the "View PDF" to get the theis with 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 sort of an issue in MIM function c2dn() implementation, which 'moves' the fc point especially at low fc values. *IIRC, same type of error happens when changing the fs.*
Analog prototype LPF is calculated as:
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:
fs=44100;
mimpim = 'mim';
samplingtime = 1/fs;
numofsamples = 2048;
[Dz1] = c2dn(Analog, samplingtime, mimpim, 1, numofsamples);
%[a, b, T] = tfdata(Dz1); %coeffs
where c2dn() is the function from MIM paper.
Plotting:
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've found out that the fc error is quite linear (I picked -3dB points for few fc and then calculated the trend f(x) (R² was 0.99989...)) so, one could probably find the source for "error" by just looking the matlab codes for MIM function c2dn().
Here's plot showing the effect of my correction f(x):
Correction formula f(x) = 0.9704384746x - 51.6991952722 works well for fc in range 60Hz-Nyqvist (fs=44.1kHz) and can even be improved a bit. (Note: correction for fc=10Hz showen in plot is done using another formula so you can omit it as a result of this f(x)).
Any suggestions regarding the source for this "error" ?
EDIT: I quess the functionality behind the issue can be located to cepstral processing block ... there's a line with remark: "% minimum phase sequence r^mn" where h is recalculated using lmn matrix created earlier in "homomorphic filtering" block:
h = h.*lmn';
When this line is turned to a comment, magnitude plots looks good for any fc (if the numofsamples parameter has suitable (big) value):
Problem is that commenting the line leads to the aliasing issue which is common for this type of filter... :( .
Any thoughts if this (homomorphic filtering/cepstral processing) can somehow be improved to result the fc point untouched?
EDIT2:
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?
