JuliaDSP / DSP.jl

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

There should be a `plot` function defined for `FilterCoefficients` #203

Closed standarddeviant closed 6 years ago

standarddeviant commented 6 years ago

I find myself writing this code and thinking, "wait, where did I write it the last time?", which suggests (to me) that it might be a good idea to have a plotting function in DSP.jl

I started writing a simple example, but kept making the plot prettier and ended up with this:

import Plots
function Plots.plot(fc::FilterCoefficients, w=linspace(0, pi, 500); fs=2, 
                    domag=true, dophase=true, 
                    doimp=false, impn=500, 
                    dostep=false, stepn=500, size=(800, 600))
    # calculate plottable vectors
    if domag || dophase
        h = freqz(fc, w)
    end
    if doimp
        ir = impz(fc, impn)
    end
    if dostep
        sr = stepz(fc, stepn)
    end

    # make vector of plots
    theplots = Plots.Plot[]
    hz = w ./ pi .* (fs/2)
    if domag
        mag = 20*log10.(abs.(h))
        ylims = [maximum(mag)-70, maximum(mag)+10] 
        push!(theplots, plot(hz, mag, lw=2, c=:blue, ylims=ylims, ylabel="dB (mag)"))
    end
    if dophase
        phase = unwrap(angle.(h))
        push!(theplots, plot(hz, phase, lw=2, c=:red, ylabel="rad (phase)"))
    end
    if doimp
        push!(theplots, plot((0:impn-1)./fs, ir, lw=2, c=:green, ylabel="imp resp"))
    end
    if dostep
        push!(theplots, plot((0:stepn-1)./fs, sr, lw=2, c=:purple, ylabel="step resp"))
    end

    # add a title
    plot!(theplots[1], title="filter response")
    # add xlabel to the last plot for space saving
    plot!(theplots[end], xlabel="Hz (freq-plots) or Seconds (time-plots)")

    plot(theplots..., layout=(sum([domag, dophase, doimp, dostep]), 1), legend=false, size=size)
end

To demo, this you can try:

using DSP
fs = 10
bpf = digitalfilter(Bandpass(1, 2; fs=fs), Butterworth(4))
plot(bpf, doimp=true, dostep=true, size=(900,500))

This feels like a good compromise between code complexity and features, but please make comments and suggestions!

ararslan commented 6 years ago

Plotting packages should be at or near the top of the package stack; that is, the vast majority of packages should not depend on a plotting package, including this one. Perhaps StatPlots.jl wants DSP functionality. You could ask there, but in the meantime we have no intention of adding the functionality here so I'm going to close this.

ssfrr commented 6 years ago

This should be doable with a plot recipe though, right? then it just has to depend on RecipesBase.jl

standarddeviant commented 6 years ago

Hi @ararslan That's a very reasonable point, and close to the response I expected.

That said, where should a FilterCoefficients aware plotting function live? I feel like visualizing filter responses is probably out of scope for StatPlots.jl, but I'll ask.

In that event that StatPlots.jl isn't a suitable home for this functionality (highly-likely), would a new repo called DSPPlots.jl or maybe DSPPlotRecipes.jl make sense?

Julia is wonderful at using custom types efficiently, but then where should visualization code live for these custom types?

standarddeviant commented 6 years ago

@ssfrr I didn't know about RecipesBase.jl, but it sounds promising based on the README.md.

So that shifts the question to: Is it acceptable to have DSP.jl depend on RecipesBase.jl? If not, creating DSPPlots.jl (or a better name without two sequential captial Ps :-) ) seems like a reasonable approach.

In the meantime, I'll read up on RecipesBase.jl.

mbaz commented 6 years ago

I don't think DSP.jl should depend on any particular plotting packages. I already use Gaston for plotting, and it'd be inefficient to bring in an entire different plotting package when I use DSP.jl. My vote would be for a DSPPlots.jl package that specializes on plotting the types and structures used by DSP.jl. Ideally, it'd be backend-agnostic.

ssfrr commented 6 years ago

The nice thing about RecipesBase is that it's super lightweight (no dependencies), the idea is to allow other packages to add plotting support without needing to depend on the actual plotting package.

I think unless it adds to the precompile time it doesn't seem like there's a downside. If it slows things down at all though I'd see the logic in separating into a DSPRecipes or DSPPlots package.

As a side note @standarddeviant DSP has a amp2db function you can use rather than 20log10 I think people who do this a lot immediately recognize the 20log10 idiom, but amp2db might be more clear. (super nitpick).

standarddeviant commented 6 years ago

That's a great point about amp2db. :relaxed:

I'll try to get the above working as a recipe and go from there.

standarddeviant commented 6 years ago

@mbaz Isn't the main point of Plots.jl to provide a back-end agnostic plotting interface?

piever commented 6 years ago

While the general principle is that defining a plot recipe only requires a tiny dependency (RecipesBase) and many packages have accepted that path (it doesn't really cost much and it gains you some functionality), if for whatever reason the package maintainers do not accept to take this dependency, the next logical thing is to add the dependency to StatPlots (that's actually the reason why StatPlots exists). In particular DSP is sufficiently lightweight that I for one am in favor of taking the dependency and would be happy to review a recipe PR to StatPlots (we should probably continue the discussion over there).

mbaz commented 6 years ago

@standarddeviant Yes, but it does not support all backends (e.g. Gaston). I also think it's too large a dependency to bring in and force all DSP users to load.

ssfrr commented 6 years ago

@mbaz as I said before I agree that DSP shouldn't depend on Plots. The question is whether it's reasonable to depend on RecipesBase, which is very lightweight.

I think the best way forward is to refactor the code above as a recipe, and then try a @time using DSP (after using some other package so we're not timing the generic using precompilation) both with and without the inclusion of RecipesBase. If the dependency doesn't meaningfully increase package load time then we can just include it right here in DSP, and if it's a burden then it should go into a separate package.

mbaz commented 6 years ago

@ssfrr I'm not opposed to depending on RecipesBase. Your proposal makes sense to me.

standarddeviant commented 6 years ago

Below is a version of the original function using plot recipes. If it looks good, I'm happy to submit a PR.

using DSP, RecipesBase
# This is all we define.  It uses a familiar signature, but strips it apart
# in order to add a custom definition to the internal method `RecipesBase.apply_recipe`
@recipe function plot(fc::FilterCoefficients; fs=2,
                    w=linspace(0, pi, 500),
                    hz=zeros(Float32, 0),
                    domag=true, dophase=true, 
                    doimp=false, impn=500, 
                    dostep=false, stepn=500, size=(800, 600))
    # markershape --> :auto        # if markershape is unset, make it :auto
    # markercolor :=  customcolor  # force markercolor to be customcolor
    grid      :=  true
    layout    :=  (sum([domag,dophase,doimp,dostep]),1)
    linewidth -->  2
    legend    :=  false

    the_xlabel = "Hz (Freq-Domain) or Seconds (Time=Domain)"

    # determine the frequencies based on w or hz
    if length(hz) > 0
        # set w according to hz
        w = hz ./ (fs/2) .* pi
    else
        hz = w ./ pi .* (fs/2)
    end

    # calculate plottable vectors
    if domag || dophase
        h = freqz(fc, w)
    end
    if doimp
        ir = impz(fc, impn)
    end
    if dostep
        sr = stepz(fc, stepn)
    end

    # make the plots
    if domag
        mag = amp2db.(abs.(h))
        @series begin
            title     :=  "filter response (fs=$fs)"
            linecolor --> :blue
            ylims     --> [maximum(mag)-70, maximum(mag)+10] 
            ylabel    --> "dB (mag)"
            if (!dophase && !doimp && !dostep)
                xlabel --> the_xlabel
            end
            hz, mag
        end
    end
    if dophase
        phase = unwrap(angle.(h))
        @series begin
            linecolor --> :red
            ylabel    --> "rad (phase)"
            if (!doimp && !dostep)
                xlabel --> the_xlabel
            end
            hz, phase
        end
    end
    if doimp
        @series begin
            linecolor --> :green
            ylabel    --> "impulse"
            if (!dostep)
                xlabel --> the_xlabel
            end
            (0:impn-1)./fs, ir
        end
    end
    if dostep
        @series begin
            linecolor --> :purple
            ylabel    --> "step"
            xlabel    --> the_xlabel
            (0:stepn-1)./fs, sr
        end
    end
end