JuliaDSP / DSP.jl

Filter design, periodograms, window functions, and other digital signal processing functionality
https://docs.juliadsp.org/stable/contents/
Other
379 stars 109 forks source link

Use different types for continuous time and discrete time filters? #84

Closed goretkin closed 3 years ago

goretkin commented 9 years ago

The documentation here appears to be using s=z^-1, which I haven't seen before. In signal processing and controls, I have always seen s to mean the Laplace-transform domain variable. Perhaps something like D=z^-1 can be used, where D is called the "delay operator".

It might also be confusing that z is used as the roots of the numerator polynomial of a transfer function and also the z-transform domain variable as in H(z)

jfsantos commented 9 years ago

We should definitely not use s=z^-1, as s has a clear "standard" meaning in signal processing.

In the case of the polynomials, z and p are used to represent the zeros and poles, respectively, but I agree that may be a bit confusing.

Would you like to send a pull request proposing a fix to these issues?

simonster commented 9 years ago

@goretkin does #97 fix your concerns?

goretkin commented 9 years ago

I am not a signal processing expert, but here are my thoughts!

I think it's still unclear (or at least I am confused) about which filters are discrete time or continuous time. The library seems to use the terms "digital" and "analog" respectively, but I think it's more accurate to describe it as DT and CT. (I see now that this is how MATLAB describes things, but I dislike it. In principle, and people do it too, you can make a sampled (discrete) time filter where signal values are represented in analog).

Even though all the implementations of filters are discrete-time filters, DSP.jl lets you design continuous-time filters and then convert them into discrete-time filters by choosing a sample rate (this is my understanding). So in terms of design, it is possible to have continuous time filters.

Whether the filter is CT or DT, an LTI, Single-Input,Single-Output filter (that starts "at rest) is completely specified by a transfer function which is a ratio of two univariate polynomials. And the different filter types (ZPK, TF, SOS) can be seen as ways of factoring or representing the transfer function, but the mathematical object they describe is the same: ratio of two polynomials.

But when Butterworth returns a ZPKFilter, the transfer function represented by the ZPKFilter is a continuous-time transfer function. Usually you'll write that transfer function as H(s)=numerator(s)/denominator(s), but it's really just a function, so you could have replaced s with something else. Now when you go to apply the filter, it's important to know what the function represents. You'll probably need to convert it to a discrete-time filter, for example, if the meaning of the transfer function is H(s) instead of H(z).

I tried delving in, but got confused by where/whether this conversion is happening. To use a Butterworth filter in filt!, it looks like the conversion goes from ZPKFilter to TFFilter to SOSFilter, but from what I understand SOS is used only for discrete-time IIR filters, so it seems like by that point the filter should have been a discrete-time filter. Furthermore, I think the documentation for BiquadFilter is what most confuses me. I think only the second definition should be included.

.. math:: H(z) = \frac{\verb!b0!+\verb!b1! z^{-1}+\verb!b2! z^{-2}}{1+\verb!a1! z^{-1} + \verb!a2! z^{-2}}

I just learned about what a Biquad filter is from wikipedia, but it seems like it is only a discrete-time filter used as the blocks of an SOS filter. If its transfer function is written down, the transform variable seems like it should definitely be z.

I am not sure if this is how the library is designed or intended, but I can see it being plausible to have, for example, ZPKFilter{:CT} and ZPKFilter{:DT}. The fields of this class will still be an array of zeros, roots, and a gain, but the zeros and roots will have different meanings depending on the parameter. And then you can use the type to let you know how to interpret the polynomial instead of having, for example, freqs and freqz. Anyway, the point of my saying this is that perhaps the least confusing thing to do in the definition of, say, ZPKFilter, is to use x as the polynomial variable. If ZPKFilter is really just a transfer function object, and that could be a transfer function for discrete-time or a continuous-time filter, then I think it would only be confusing to use one of s or z.

simonster commented 9 years ago

"Digital filter" and "analog filter" seem to be the terms used by most of the references I've seen. I guess "discrete-time filter" is a broader term; Wikipedia notes that there are analog discrete-time filters that use switched capacitors.

It is true that all of the filter representations are different factorizations of the transfer function, and so calling the polynomial representation the transfer function isn't quite right, but this is the terminology used by MATLAB and others.

At some point I changed the filter code to avoid the intermediate TF conversion during filter design. Now all of the design happens in ZPK and the ZPK->SOS conversion happens only when filt is called.

I know next to nothing about analog filters, but a quick Google shows a lot of references that imply that analog SOS filters are commonly used because they are less sensitive to component variation.

A biquad is just a TFFilter with 3 coefficients. It's a separate type from TFFilter only because that allows more efficient filter implementation.

We support analog/continuous time filter design only because it's simple to implement, but I don't think it's the primary focus of this library. I guess having a separate continuous time filter object would stop someone from trying to implement an analog filter on discretized data and would let us use multiple dispatch for freqs/freqz. There are also some subtleties to how we should interpret the coefficients of CT vs. DT TFFilters if the two coefficient vectors are different lengths that we don't currently capture.

I have no objections to using x for ZPKFilter, although it looks like MATLAB only uses s in most of its docs.

simonster commented 9 years ago

Updated the title to reflect what I think is the only remaining question here.

kumarbal commented 4 years ago

I have a related issue I just opened up in https://github.com/JuliaDSP/DSP.jl/issues/341. I think it is great that you support continuous-time filtering, but I would prefer that the connection to sampled data representations be handled outside of the basic design package. There is something broken with the interface to filter as I expect the Laplace complex frequency to be convertible to j\omega. Also, I don't understand the scaling of the filter response when using the freqs function and wuld like independence from the sampling frequency.

An analog filter should be convertible to a sampled data representation or a digital filter.

martinholters commented 3 years ago

Fixed by #284.