KVR Audio

juha_p · Post by **juha_p** » Sat Jul 30, 2016 8:14 am

How to plot the difference (error) between analog and digital filter in Octave/Matlab?

Example code:

fs = 44100;

%% analog model
a2 = [0.000015 1];
b2 = [0.000050 1];
Ds = tf(a2, b2);
Ds = Ds/dcgain(Ds);

%% digial 1
B1 = [1 -0.62786881719628800 0];
A1 = [-0.08782333709141920 0.45995451989513100 0];
Dz = tf(A1, B1, 1/fs);
Dz = Dz/dcgain(Dz);

%% digital 2
B2 = [1 0.8868033819284362 -0.312597640666886 -0.3462305068062579 -0.04427163895789011];
A2 = [0.4601686131709282 0.6087942052682088 0.1563216019782706 -0.03386380454131521 -0.007717020378690522]; 
Dz2 = tf(A2, B2, 1/fs);
Dz2 = Dz2/dcgain(Dz2);

wmin = 2 * pi * 20; % 20Hz
wmax = 2 * pi * (fs/2); 

bode(Ds, 'r', Dz, 'k', Dz2, 'g', {wmin, wmax});
legend('Ds', 'Dz1', 'Dz2','location', 'northwest');

% plot errors

I tried something but got errors related to different sample rate (not given for analog model) or vector size (JOS's method).

Miles1981 · Post by **Miles1981** » Sat Jul 30, 2016 8:29 am

I have that in Python for my library. For instance, you use scipy.signal.freqs to compute the analog continuous spectrum and scipy.signal.freqz for the digital one (https://github.com/mbrucher/AudioTK/blo ... _ATK_LP.py for a complete example of freqz with custom sizes).

juha_p · Post by **juha_p** » Sat Jul 30, 2016 12:04 pm

Hmm...

Code: Select all

fs = 44100;
n = 512;

a2 = [0.000015 1];
b2 = [0.000050 1];

w = logspace(0,4,n);
Ha = freqs(b2, a2, w);

B1 = [1 -0.62786881719628800 0];
A1 = [-0.08782333709141920 0.45995451989513100 0];

[H0 w] = freqz(B1, A1, n);
[Bh, Ah, Sig0] = invfreqz(H0, w, 2, 2);
[Hh, wh] = freqz(Bh, Ah, n);

plot(w,[abs(Ha) abs(Hh)])
xlabel("Fs (rad/sample)");
ylabel("magnitude");
legend('analog','digital');

results an error in plot() command:

Code: Select all

error: horizontal dimensions mismatch (1x128 vs 128x1)

Any thoughts?

stratum · Post by **stratum** » Sat Jul 30, 2016 1:05 pm

error: horizontal dimensions mismatch (1x128 vs 128x1)

transpose?

juha_p · Post by **juha_p** » Sat Jul 30, 2016 5:00 pm

Yes, transposing the analog filter clears that error.

@Miles1981
There's a this line plt.title('Digital filter frequency response') in your linked source code ... does it mean that you just plot the various curves to a same plot ? If it does than I do it already in my linked code. If you plot the errors agains analog model there then do you mean I would need to use functions loglog() and semilogx() in octave code?

I'm after some instructions for plotting the error between the analog model and the digital filter representing it. This is so easy to do in Excel/Calc ... is it really this complicated in Octave (Matlab)?

juha_p · Post by **juha_p** » Sat Jul 30, 2016 6:46 pm

OK, found one possible solution for to plot magnitude errors. Solution is based in code listing given in 1st post (see how the gain 'normalization' is done there).

Miles1981 · Post by **Miles1981** » Sat Jul 30, 2016 6:48 pm

I just plotted some magnitude results on the same graph to compare different EQs. But you can do different things.
Once you have the results from freqs and freqz, you cqn just substract and square the result to get the error.

juha_p · Post by **juha_p** » Sat Jul 30, 2016 7:00 pm

Miles1981 wrote:I just plotted some magnitude results on the same graph to compare different EQs. But you can do different things.
Once you have the results from freqs and freqz, you cqn just substract and square the result to get the error.

Thanks. I'll try that way aswell.

How to plot the difference between analog and digital filter in Octave/Matlab