DSP.jl how to use ZeroPoleGain with specifications in Hz?

General Usage
1

Let’s say I want to implement a digital filter that is defined by its zeros and poles to filter signals sampled at fs Hz. It seems logical to use DSP.Filters.ZeroPoleGain(z,p,k). But there is little information about the units for z and p. If the zeros and poles are in Hz+imHz, how do I obtain suitable values for z and p? If zHz is a pole in Hz+imHz, do I have to divide it by fs to obtain z? Should I multiply by 2*pi?. Map from the s-plane to the z-plane? Sorry for the stupid question.

1 Like
2

You may need to convert the frequencies in Hz to rad/sample.
Pi rad/sample equals to the Nyquist frequency ( = sampling frequency / 2 = fs/2 in Hz).
Then the frequency vector frs in rad/sample which we need to feed into freqresp(filter, frs) is given by frs = pi * fhz / (fs/2), where fhz are the frequencies in Hertz.

3

Thanks, but it does not work as expected. Here is my code:

using DSP, Plots
#4 zeros in Hz
zhz=[0.,0.,0.,0.]
#8 poles in Hz (2 of magnitude 1Hz and 6 of magnitude 20 Hz)
phz=[-0.707(1-im), -0.707(1+im), -19.27-5.16im, -19.27+5.16im, -14.11(1-im), -14.11(1+im), -5.16-19.27im, -5.16+19.27im]
fs=44100. #sampling frequency
f=ZeroPoleGain(2π/fs*zhz,2π/fs*phz,1)
fhz=0.01:0.01:400.
frs=2π/fs*fhz
A=20log10.(abs.(freqresp(f,frs)))
plot(fhz,A)

This filter should feature a minimum attenuation at 10 Hz and increasingly attenuate as the distance to 10 Hz increases. Here the attenuation I get is more or less constant.

4

Sorry, I don’t understand what you mean by the “20 Hz magnitude” poles. Some of the poles are outside the unit circle and the filter will not be stable, having exponentially increasing components.

I recommend that you try to reproduce using DSP.jl a well documented example, like this one.

5

After multiplication by \pi/(fs/2) my poles and zeros originally specified in Hz will be inside the unit circle.

The example you suggested works.

using DSP, Plots
z=[0.95+0.31im, 0.95-0.31im]
p=[0.9-0.3im, 0.9+0.3im]
f=ZeroPoleGain(z,p,1)
frs=0:0.001:π
A=20log10.(abs.(freqresp(f,frs)))
plot(frs,A)

But that does not tell me what to do with my poles and zeros in Hz.

6

I don’t think it is valid to apply frequency conversions directly to the poles and zeros. They should be applied only to the frequency vector input to the freqresp() function, that is expected in rad/sample for a digital filter.

NB:
and then plot the filter frequency response in a Hertz scale with plot(fhz, A)

8

In that example, where there is one pole at: p1 = 0.9 + 0.3*im, the “corresponding frequency” (in the transfer function) in Hertz can be computed using:

angle(p1)/π * fs/2
9

Thanks. Then DSP.Filters.ZeroPoleGain() takes arguments in the z-plane. And I have to map my zeros and poles in the s-plane onto the z-plane. A possible approach is to use the bilinear transform. Do you know of any Julia equivalent to the bilinear() Matlab function that is used here:

p=[ -0.707+j* 0.707
    -0.707-j* 0.707
   -19.27+j* 5.16
   -19.27-j* 5.16
   -14.11+j*14.11
   -14.11-j*14.11
    -5.16+j*19.27
    -5.16-j*19.27];
%zeros
z=[0 0 0 0]';
%scalar gain
k=10^(116/20);
%get analog filter coefficients
[bG,aG] = zp2tf(z,p,k);
%get filter response
[hG,fG] = freqs(bG,aG,[0 logspace(log10(0.2),log10(200),1000)]);

%%Try to get digital coefficients from bilinear transform
[bz,az]=bilinear(bG,aG,fs);
%get digital filter response
[Hz, fz]= freqz(bz,az,100*fs,fs);
```

From https://se.mathworks.com/matlabcentral/answers/6373-s-plane-to-z-plane-transformation.
11

There is a DSP.Filters.bilinear() function that is not documented. See this other post and the github source code.