Tweaked One pole LPF
![Post Post](./styles/Bl%40ck/imageset/icon_post_target.gif)
FYI, not long ago I got an idea of trying to speed up MZT/IIM based one pole LPF by using Taylor series to approximate exp() function. After some tryouts I found out that use of just degree 2 polynomial resulted approximation error which actually improved the resposne of the MZT/IIM filters when cut-off frequency was approaching Nyqvist. Not perfect filter but better ... and probably can be improved more but my math skills are not that good... .
So, as winked, coefficients are calculated this way (Matlab/Octave code):
Comparison against MZT, IIM and BLT: https://i.postimg.cc/zDPyPrRz/lpf-robo.png (cut-offs: 1k, 5k, 10k, 15k, 20k)
I also did run speed comparison against the original exp() based filter implementation and got about 20x speed up with this tweak (Ubuntu, g++, --ffast-math enabled).
Octave code for the linked plot image:
So, as winked, coefficients are calculated this way (Matlab/Octave code):
- Code: Select all
// fs - in Hz (tested only @ 44100)
// fct - in µs, if Hz used then convert to µs with formula 1/(2*pi*fc)
x = 1.0/(fs*fct);
a0 = 1.0;
a1 = -(1 + x + 0.5 * x^2); // Taylor
b0 = a0 + a1;
b1 = 0.0;
a = [a0 a1];
b = [b0 b1];
Comparison against MZT, IIM and BLT: https://i.postimg.cc/zDPyPrRz/lpf-robo.png (cut-offs: 1k, 5k, 10k, 15k, 20k)
I also did run speed comparison against the original exp() based filter implementation and got about 20x speed up with this tweak (Ubuntu, g++, --ffast-math enabled).
Octave code for the linked plot image:
- Code: Select all
% --------------------
% OCTAVE PACKAGES
% --------------------
pkg load signal
% --------------------
%clf
fs=44100;
fc=1000;
fct=1/(2*pi*fc);
fs2 = fs/2;
N=1; % order
% -----------------------------
% Analog model
% --------------------
w0 = 2*pi*fc;
Analogb = 1;
Analoga = [1 w0];
Analog = tf(Analogb, Analoga);
Analog = Analog/dcgain(Analog);
% --------------------
% MZT
% --------------------
p1 = exp(-1.0/(fs*(1/(2 * pi * fc))))
z1 = 0.0;
MZTa0 = 1.0;
MZTa1 = -p1;
MZTb0 = 1.0-p1;
MZTb1 = -z1;
MZTa = [MZTa0 MZTa1];
MZTb = [MZTb0 MZTb1];
MZT=tf(MZTb, MZTa, 1/fs);
% --------------------
% IIM
% --------------------
[IIMb,IIMa]=impinvar(Analogb,Analoga,fs);
IIM=tf(IIMb, IIMa, 1/fs);
IIM=IIM/dcgain(IIM);
% --------------------
% MZT approximated
% --------------------
x = 1.0/(fs*fct);
MZTa0 = 1.0;
MZTa1 = -(1 + x + 0.5 * x^2);
MZTb0 = MZTa0 + MZTa1;
MZTb1 = 0.0;
MZTa = [MZTa0 MZTa1];
MZTb = [MZTb0 MZTb1];
MZT2=tf(MZTb, MZTa, 1/fs);
% --------------------
% BLT
% --------------------
w0 = 2*pi*(fc/fs);
%BLTb0 = sin(w0);
%BLTb1 = sin(w0);
%BLTa0 = cos(w0) + sin(w0) + 1.0;
%BLTa1 = sin(w0) - cos(w0) - 1.0;
%
%BLTa = [BLTa0 BLTa1];
%BLTb = [BLTb0 BLTb1];
%
%BLT=tf(BLTb, BLTa, 1/fs);
%BLT=BLT/dcgain(BLT);
BLT = c2d(Analog, 1/fs, 'prewarp', w0);
% --------------------
% Plot
% --------------------
nf = logspace(0, 5, fs2);
figure(1);
% analog model
[mag, pha] = bode(Analog,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'g', 'linewidth', 2);
axis([10 fs2 -30 1]);
hold on;
%semilogx(nf, pha, 'color', 'g', 'linewidth', 2, 'linestyle', '--');
% BLT
[mag, pha] = bode(BLT,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'm', 'linewidth', 1.0, 'linestyle', '-');
%semilogx(nf, pha, 'color', 'm');
% MZT
[mag, pha] = bode(MZT,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'k', 'linewidth', 1.0, 'linestyle', '-');
%semilogx(nf, pha, 'color', 'k', 'linewidth', 1.0, 'linestyle', '--');
% IIM
[mag, pha] = bode(IIM,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'c', 'linewidth', 1.0, 'linestyle', '--');
%semilogx(nf, pha, 'color', 'c', 'linewidth', 1.0, 'linestyle', '--');
% MZT approximated
[mag, pha] = bode(MZT2,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'r', 'linewidth', 2.0, 'linestyle', '--');
%semilogx(nf, -pha, 'color', 'r', 'linewidth', 2.0, 'linestyle', '--');
grid on;
str=num2str(fs);
str2=num2str(fc);
str3=num2str(N);
str = sprintf("LPF (various TF), order:%s, fs=%s, fc=%s, ",str3, str, str2);
title(str);
legend('Analog', 'BLT', 'MZT', 'IIM', 'MZTapprox', 'location', 'southwest');
xlabel('Hz');ylabel('dB');