Octave/Matlab freqz in C++

[rd2play]

I am learning DSP on my own and it's slow going, but it's fun. I wish I went for it in college... It must be great to have a question and get the answer next class, instead of going through online resources and books and piece together the answer. You could say you learn a lot this way, but it sure is frustratingly slow at times.

Anyway, I have a basic LP DFI filter in Octave that plots nicely the frequency response. I put the code in C++ and I get the same output y(n) (obviously). However, this is when I don't know what to do next.

In Octave I simply do:

Code: Select all

[h,w] = freqz(y,1,4096,Fs);

and then I do:

Code: Select all

semilogx(w,20*log10(abs(h)));

The question is what do I do in C++ for "freqz"? What do I do with "y"?

I have this code, but not sure how to get to a plot... (I already calculated the coefficients earlier to get y(n))

Code: Select all

double num = (b0 * b0) + (b1 * b1) + (b2 * b2) + (2.0 * cosW0 * ((b0 * b1) + (b1 * b2))) + (2.0 * cos2W0 * b0 * a2);
double den = 1.0 + (a1 * a1) + (a2 * a2) + (2.0 * cosW0 * (a1 + (a1 * a2))) + (2.0 * cos2W0 * a2);

BTW, I'm using Dislin for plotting. I kinda got a hang of it, but still working on it

[mystran]

If you have some transfer function (b0+b1*z^-1+b2*z^-2)/(1+a1*z^-1+a2*z^-2) then this is simply the z-transform (similar to the Laplace transform of continuous time systems) and the discrete-time Fourier transform is found on the unit-circle of the z-transform z-plane (similar to how Fourier transform is found on the Laplace imaginary axis).

So in order to evaluate the Fourier transform (ie. frequency response) at some frequency w=2*pi*f/fs, all we have to do is take z=exp(i*w)=cos(w)+i*sin(w) to get the relevant point on the unit circle and compute the value of the transfer function at this point. If you then transform the resulting complex frequency response into polar form, the complex magnitude gives you the amplitude response (which you can then convert to logarithmic scale if desired) and the complex argument gives you the phase response.

In C++ you can simply do something like this:

Code: Select all

typedef std::complex<double> complex;
double pi = acos(-1.);
complex zm1 = std::exp(complex(0, -2*pi*f/fs));
complex zm2 = zm1*zm1;
double magnitude = std::abs((b0+b1*zm1+b2*zm2)/(1.+a1*zm1+a2*zm2));

A real-only expression for computing just the magnitude response can be obtained by expanding this complex expression in terms of real variables and simplifying, but that expression will only ever make any sense once you first understand the complex version.

[rd2play]

Thanks, I appreciate the help. I understand the theory and formulas of Laplace, Fourier, z-transform, but going from that theoretical understanding to practical applications is where I miss a formal DSP education or at least access to some sort of tutoring.

I'm trying to understand the code... Why zm1 * zm1?

freqz uses the output y to calculate the response. Where does y fit in?

[mystran]

I'm trying to understand the code... Why zm1 * zm1?

freqz uses the output y to calculate the response. Where does y fit in?

If you look carefully at the argument for std::exp, there's a negative sign. What I did there was merge the exponent of z^-1 with the 2*pi*f/fs term, which is valid because (x^a)^b=x^(a*b) holds for complex numbers. Since z^-1 is not a valid C++ identifier, I called it zm1 instead (ie. "minus 1").

The reason for doing this is that then we can work with positive powers of z^-1 (=1/z) rather than negative powers of z. Since (x^a)*(x^b)=x^(a+b) also holds, zm1*zm1 is simply (z^-1)*(z^-1)=z^-2 which is what zm2 stands for.

The freqz function computes the filter's response from the transfer function coefficients. If we take a digital filter and feed it with a unit impulse, then for a FIR filter (ie. I think that's what passing 1 as the second vector does) the output is simply the original coefficient vector. For an IIR filter the response is infinite (hence the name) so we'll have to truncate and if we treat a finite length recording as a FIR filter then this will approximate the response of the original IIR with the approximation improving as we make the FIR longer (there are better ways to approximate an IIR with a FIR though).

I'm not really an Octave/Matlab user, but I'd imagine that if you give freqz some recording of the output of the filter, then freqz will likely just treat it as a coefficient vector of a FIR and if this FIR is long enough then it will look fairly similar to the response of whatever IIR you might have started with. You could also give it your original a,b coefficients and then it would basically compute what I did above.

[rd2play]

Since (x^a)*(x^b)=x^(a+b) also holds, zm1*zm1 is simply (z^-1)*(z^-1)=z^-2 which is what zm2 stands for.

Now that you explained it, all I can say is duh... I've seen the formula countless times, but never in computer language. I didn't even know until now that there is a complex data type. (This led me right away to the question "how does C deal with complex numbers?" and I found the answers online; interesting... I'll have to dig a little deeper)

I agree with your explanation/assumption of freqz. I kinda assumed something similar, but I just don't trust myself. I have to find an online tutor. Thank you!

[mystran]

Since (x^a)*(x^b)=x^(a+b) also holds, zm1*zm1 is simply (z^-1)*(z^-1)=z^-2 which is what zm2 stands for.

Now that you explained it, all I can say is duh... I've seen the formula countless times, but never in computer language.

Another way to explain the same thing is that if z^-1 stands for 1/z then (1/z)*(1/z)=1/(z^2)=z^-2.

Personally, in code I often actually just write z when I mean z^-1 and z2 when I mean z^-2, because really most of the time you're just treating z^-1 as the main variable and you don't need "z" as a variable for anything, but here I wanted to try to make things a bit less confusing for pedagogical reasons.

[rd2play]

Got it. I wasn't confused about the notation itself, just didn't understand why they were multiplied. Obviously to get z^-2...

[rd2play]

For the sake of completeness, I found this interesting bit of information. The code in blue below returns the same frequency response as Matlab's freqz (the code in red). Pretty cool...

The code in blue is adapted from the book "Digital Signal Processing, Principles and Applications" by Thomas Holton (page 157).

===========================================
x = [1;zeros(255,1)];
Fs = 48000;
f = 1000;
Q = 0.707;
dBGain = 1;
N = 4096;

y = somebiquadfunction(x,Fs,f,Q,dBGain);
[hh,ww] = freqz(y,1,N,Fs);

b = y ;
a = 1;
lenb = (0:length(b) - 1);
lena = (0:length(a) - 1);
w = pi * linspace(0,1 - 1/N ,N) ' ;
h = (exp(-1j*w*lenb)*b(:)) ./ (exp(-1j*w*lena)*a(:)) ;

figure
subplot(2,1,1)
plot( w / pi , abs ( h ))
subplot(2,1,2)
plot(ww / pi , abs ( hh ))

[mystran]

b = y;
a = 1;
lenb = (0:length(b) - 1);
lena = (0:length(a) - 1);
w = pi * linspace(0,1 - 1/N ,N) ' ;
h = (exp(-1j*w*lenb)*b(:)) ./ (exp(-1j*w*lena)*a(:)) ;

Let's look at this code very carefully and try to see what is really going on.

b here is the vector of numerator coefficients, while a is the vector of denominator coefficients. Then lenb and lena vectors of sequential indexes (ie. [0,1,2,3..], with the same length as b and a respectively.

linspace generates a vector w of linearly spaced (normalized) frequencies for us; these are the frequencies we are evaluating, though we could also explicitly choose a set of specific frequencies.

Finally it computes the frequency response for each of the frequencies w basically as:
(b[0] + b[1]*exp(-i*w*1) + b[2]*exp(-i*w*2) ... )
/ (a[0] + a[1]*exp(-i*w*1) + a[2]*exp(-i*w*2) ...)

MATLAB people love abusing their operators to try and write code that makes it non-obvious what is going on, but that's what it's really computing.

If that makes sense, then we observe that you're not really giving this function the coefficients of the biquad (ie. b=[b0,b1,b2] and a=[1,a1,a2]), but rather sampling a truncated impulse response y and then pretending that these are the numerator coefficients of a FIR (ie. denominator a = [1, 0, 0, ...]).

At higher frequencies this FIR will approximate the original IIR biquad (it'll approximate it everywhere, but as we go down in frequency the approximation will be increasingly poor) and taking a longer impulse response before truncation can improve the approximation, but we're really computing the frequency response of this truncated IR, even though both freqz and the code in blue are designed to compute the frequency response of an IIR directly.

Computing the frequency response of a FIR also works, because a FIR is just a special case of an IIR and this approach (well, sort of; usually you'd use an FFT which is faster to compute) is what you use to measure filters and systems when you don't know the coefficients... but if you have a biquad and you know the coefficients, then with b=[b0,b1,b2] and a=[1,a1,a2] you will get the same results as the C++ code example I posted above (which could also be expanded to compute frequency responses for a FIR just by computing as many numerator terms as you have and dropping the denominator).

Does this distinction make sense now?

[rd2play]

Thank you, it totally makes sense. You had mentioned this before... I didn't know that you can use the IIR output as an FIR numerator, to approximate the IIR.

Funny about the Matlab operators... I found out the hard way that a "." or a " ' " make a huge difference... I just got started with Matlab, so I'm sure there is a lot I don't know, but those are my attempts to plot the 3 options. It looks like the second one, the one using coefficients in the blue code, has a slight gain. The mouse pointer tells me it's around 15.

Output.jpg

[mystran]

Thank you, it totally makes sense. You had mentioned this before... I didn't know that you can use the IIR output as an FIR numerator, to approximate the IIR.

Sure. Just keep in mind that since the error essentially depends on the remaining amplitude of the decaying tail around the truncation point, so higher Q (decays slower per cycle) or lower frequencies (individual cycles are longer) will make the error (potentially much) worse for a given FIR length.

I never bothered to learn MATLAB, so can't really help you with that stuff... but as a general tip when plotting response it often makes sense to set an explicit lower limit (eg. -100dB or whatever), because if you rely on auto-scaling, then there's the possibility that one of the sampling points is inside a deep notch and then you get several hundred dB worth of vertical scale which is mostly irrelevant because in practice anything below -120dB or so is essentially pure silence and the most interesting stuff to look for is usually within the first 60dB or so.

[rd2play]

Cool stuff, I really appreciate the education! What are you using if not Matlab, is it Python?

[vasyan]

b = y;
a = 1;
lenb = (0:length(b) - 1);
lena = (0:length(a) - 1);
w = pi * linspace(0,1 - 1/N ,N) ' ;
h = (exp(-1j*w*lenb)*b(:)) ./ (exp(-1j*w*lena)*a(:)) ;

Let's look at this code very carefully and try to see what is really going on.

b here is the vector of numerator coefficients, while a is the vector of denominator coefficients. Then lenb and lena vectors of sequential indexes (ie. [0,1,2,3..], with the same length as b and a respectively.

linspace generates a vector w of linearly spaced (normalized) frequencies for us; these are the frequencies we are evaluating, though we could also explicitly choose a set of specific frequencies.

Finally it computes the frequency response for each of the frequencies w basically as:
(b[0] + b[1]*exp(-i*w*1) + b[2]*exp(-i*w*2) ... )
/ (a[0] + a[1]*exp(-i*w*1) + a[2]*exp(-i*w*2) ...)

MATLAB people love abusing their operators to try and write code that makes it non-obvious what is going on, but that's what it's really computing.

Dear colleagues, please help me figure out how to translate this formula from C++ code. I have been breaking my brain for a whole week, I can not find a single implementation of this algorithm in C ++, as well as a single converter from matlab to C++.
After all, you can’t put the code on matlab into the plugin.
Sorry for my English.

[juha_p]

Could this help you out (though, its for biquad filter only)? : https://www.earlevel.com/main/2016/12/0 ... ment-53251

There's also this IT++ package which implements freqz by using fft - https://sourceforge.net/p/itpp/git/ci/m ... design.cpp

[vasyan]

Could this help you out (though, its for biquad filter only)? :

I need to get the frequency response from the convolution pulse (IR) in order to graph it, so precision is not required.
It seems that the last link is what is required. I will sort this out. Thank you