Issue with freqz() with a single complex zero/pole

[adriano-vilela]

Hello,

I'm experiencing a problem with the freqz() function when computing the frequency response of a discrete-time filter with a single complex zero (as opposed to a complex conjugate pair of zeros). For example, imagine I want to compute the frequency response of the filter H(z) = 1 - c*z^{-1}, where c is the only zero of the filter (of course, in this case, there's also a single pole at z=0). If I make c a complex number, I get a frequency response that is mirrored (i.e., flipped horizontally). The code below illustrates the problem:

using DSP, PyPlot
c = 0.5+0.5im;
num = [1, -c];
den = [1];
w = linspace(0, 2*pi, 500);
H = freqz(PolynomialRatio(num,den),w);
plot(w,20*log10(abs(H)));

The attached figures show the problem (one figure shows the output I get in Julia with DSP.jl, whereas the other figure shows the output computed with Matlab).

Does anybody know what's going on? Am I missing something here?

Thank you,

Adriano

PS: The same problem occurs with a single complex pole.

[ma-laforge]

As commented in https://github.com/JuliaDSP/DSP.jl/issues/135, your issue is actually that:

num = [1, -c];

corresponds to:

- (0.5 + 0.5im) + 1.0⋅z^-1

(the coefficients are reversed)

This can be confirmed by show()-ing the value of the numerator (which is the .b field). The following example shows how to do this (and corrects the numerator polynomial):

c = 0.5+0.5im;
num = [-c, 1];
den = [1];
p = PolynomialRatio(num,den)
@show p.b #=>p.b = Poly(1.0 - (0.5 + 0.5im)⋅x)

I am not saying this choice for polynomial ordering is either good or bad. I am simply showing you how to examine the Julia data structures created/used by DSP.jl.

[adriano-vilela]

Thanks for your reply.

We have long figured the problem was the reversed coefficients. That's why I asked my student to report the problem here (which he did in issue #135).

According to the documentation for the PolynomialRatio() function (http://dspjl.readthedocs.io/en/latest/filters.html#PolynomialRatio) the polynomial H(z) = 1 - c*z^{-1} must be entered as [1, -c] (i.e., as increasing powers of z^{-1}). This is in line with Matlab's implementation and most textbooks on Digital Signal Processing.

Thus, there's a conflict between the documentation and the implementation here. Either one must be changed in order to address the issue. Since the documentation is in line with what's most used out there, my opinion is that the implementation should be changed.

Adriano

[adriano-vilela]

Since this bug has been around for a while now, I decided to take a look at the source code to see if I could fix it myself. In doing so, I found out there are two things wrong:

1. Line 9 of the file src/Filters/response.jl shows:

ejw = exp(-im * w)

However, the correct is

ejw = exp(im * w)

This is because the frequency response is obtained simply by evaluating the system function over the unit circle by making z = e^(jω). Just correcting this solves the problem reported here. However, this also introduces some problems. For example, the code below:

num = [1, 2, 3];
den = [2, 1];
F = PolynomialRatio(num,den);

produces the following error:

ERROR: InexactError()
Stacktrace:
 [1] convert(::Type{Int64}, ::Float64) at ./float.jl:679
 [2] copy!(::IndexLinear, ::Array{Int64,1}, ::IndexLinear, ::Array{Float64,1}) at ./abstractarray.jl:648
 [3] DSP.Filters.PolynomialRatio{Int64}(::Polynomials.Poly{Int64}, ::Polynomials.Poly{Int64}) at /home/adriano/.julia/v0.6/DSP/src/Filters/coefficients.jl:41
 [4] DSP.Filters.PolynomialRatio(::Array{Int64,1}, ::Array{Int64,1}) at /home/adriano/.julia/v0.6/DSP/src/Filters/coefficients.jl:52

2. The second problem is more fundamental. Basically, using a single PolynomialRatio representation for both continuous-time and discrete-time systems leads to an ambiguity. For example, consider we construct the PolynomialRatio object below:

num = [1, 2, 3];
den = [2, 1];
F = PolynomialRatio(num,den);

In the s-domain (continuous time), this represents the following system function:

H(s) = (1⋅s^2 + 2⋅s + 3)/(2⋅s + 1) = (3 + 2⋅s + 1⋅s^2 )/(1 + 2⋅s)

Since Polynomials.Poly() expects the coefficients to be in ascending order, we can represent H(s) as a Polynomials.Poly object by simply reversing the order of the numerator and denominator:

F = Poly([3, 2, 1], [1, 2])

That's exactly what the PolynomialRatio() constructor does (see src/Filters/coefficients.jl, line 52).

Consider again the example above:

num = [1, 2, 3];
den = [2, 1];
F = PolynomialRatio(num,den);

In the z-domain (discrete time), this represents the following system function:

H(z) = (1 + 2⋅z^(-1) + 3⋅z^(-2))/(2 + 1⋅z^(-1)) = (3 + 2⋅z + 1⋅z^2)/(z + 2⋅z^2)

This can be represented as Polynomials.Poly object by doing:

F = Poly([3, 2, 1], [0, 1, 2])

However, the polynomial that gets constructed by PolynomialRatio() is

H(z) = (3 + 2⋅z + 1⋅z^2 )/(1 + 2⋅z)

which is obviously wrong.

Thus, as far as I understand, PolynomialRatio representations of discrete-time systems in DSP.jl are wrong. I don't think it's possible to use a single representation for both continuous-time and discrete-time systems. This is because continuous systems are represented as a ratio of polynomials in descending order of s, whereas discrete systems are represented as a ratio of polynomials in ascending order of z^(-1). One possible way to fix this would be to pass an extra parameter to the constructor in order to tell it whether we want a continuous-time or a discrete-time system (it could be just "s" or "z"; see the function tf() in Matlab as an example). This would allow us to compute the correct vectors (which cannot be computed by simply reversing them) that should be passed to Polynomials.Poly() in the case of discrete-time systems.

[martinholters]

Fixed by #284.