GitHub

Convolution matrix

This example illustrates 1D signal convolution represented as matrix operations (for various boundary conditions) using the Julia language.

This page comes from a single Julia file: conv-mat.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: conv-mat.ipynb, or open it in binder here: conv-mat.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "ColorSchemes"
        "DSP"
        "FFTW"
        "FFTViews"
        "InteractiveUtils"
        "LaTeXStrings"
        "LinearMapsAA"
        "MIRTjim"
        "Plots"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using ColorSchemes
using DSP: conv
using FFTW: fft, ifft
using FFTViews: FFTView
using InteractiveUtils: versioninfo
using LaTeXStrings
using LinearMapsAA
using MIRTjim: jim
using Plots: RGB, default, gui, heatmap, plot, plot!, savefig, scatter!, @layout
using Plots.PlotMeasures: px
default(markerstrokecolor=:auto, markersize=11, linewidth=5, label="",
 tickfontsize = 11, labelfontsize = 15, titlefontsize=18)

Filter

Define a 5-tap finite impulse response (FIR) filter $h[k]$ with support -2:2 and display it.

psf = [1, 3, 5, 4, 2] # filter impulse response
K = length(psf)
hp = Int((K-1)/2) # 2 (half width)

cols = [ColorSchemes.viridis[x] for x in range(0,1,K)]

pp = plot(widen=true, tickfontsize = 14, labelfontsize = 20,
    xaxis = (L"k", (-1,1) .* 4, -4:4),
    yaxis = (L"h[k]", (-0.2, 5.2)),
    size = (600, 250), left_margin = 15px, bottom_margin = 18px,
)
plot!([-5, 5], [0, 0], color=:black, linewidth=2)
for k in 1:K
    c = psf[k]
    c = cols[c]
    plot!([k-hp-1], [psf[k]], line=:stem, marker=:circle, color=c)
end
for k in (hp+1):4
    scatter!([k], [0], color=:grey)
    scatter!([-k], [0], color=:grey)
end
gui()
# savefig(pp, "psf.pdf")

Convolution matrices

Define convolution matrices for various end conditions. The names of the conditions below match those of Matlab's conv function. Submit a github issue to request elaboration on these cases.

N = 10 # input signal length
kernel = FFTView(zeros(Int, N))
kernel[-hp:hp] = psf;

'full'

for zero end conditions where the output signal length is $M = N + K - 1$.

Az = LinearMapAA(x -> conv(x, psf), (N+K-1, N))
Az = Matrix(Az)
Az = round.(Int, Az) # because `conv` apparently uses `fft`
size(Az)
(14, 10)

'circ'

for periodic end conditions for which $M = N$.

Ac = LinearMapAA(x -> ifft(fft(x) .* fft(kernel)), (N,N) ; T = ComplexF64)
Ac = Matrix(Ac)
Ac = round.(Int, real(Ac))
size(Ac)
(10, 10)

'same'

for matching input and output signal lengths, so $M = N$.

As = Az[hp .+ (1:N),:]
size(As)
(10, 10)

'valid'

where output is defined only for samples where the shifted impulse overlaps the input signal, for which $M = N-K+1$.

Av = Az[(K-1) .+ (1:(N-(K-1))),:]
size(Av)
(6, 10)

Display convolution matrices

using colors that match the 1D impulse response plot.

cmap = [RGB{Float64}(1.,1.,1.) .* 0.8; cols];

jy = (x,t,y) -> jim(x'; color=cmap, ylims=(0.5,14.5), ylabel=y, labelfontsize=12);

pz = jy(Az, "'full'", L"M=N+K-1")
pv = jy(Av, "'valid'", L"M=N-K+1")
ps = jy(As, "'same'", L"M=N")
pc = jy(Ac, "circulant", L"M=N")
# plot(pz, pv, ps, pc; size=(400,400))
p4 = plot(pz, pv, ps, pc; size=(1000,200), layout=(1,4), left_margin = 20px)

# savefig(p4, "match4.pdf")
# savefig(pz, "a-full.pdf")
# savefig(pv, "a-valid.pdf")
# savefig(ps, "a-same.pdf")
# savefig(pc, "a-circ.pdf")
Example block output
l = @layout [
 a{0.3w} b{0.3w}
 c{0.3w} d{0.3w}
]
p22 = plot(pz, pv, ps, pc) #, layout = l, size = (500,500))
# savefig(p22, "all4.pdf")
Example block output

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
11-element Vector{SubString{String}}:
 "Julia Version 1.11.1"
 "Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [aaaa29a8] Clustering v0.15.7
  [35d6a980] ColorSchemes v3.27.1
⌅ [3da002f7] ColorTypes v0.11.5
⌃ [c3611d14] ColorVectorSpace v0.10.0
  [717857b8] DSP v0.7.10
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.7.0
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.8.0
  [587475ba] Flux v0.14.25
  [a09fc81d] ImageCore v0.10.4
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.0.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.20.1
  [7035ae7a] MIRT v0.18.2
  [170b2178] MIRTjim v0.25.0
  [eb30cadb] MLDatasets v0.7.18
  [efe261a4] NFFT v0.13.5
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.5
  [429524aa] Optim v1.10.0
  [91a5bcdd] Plots v1.40.9
  [f27b6e38] Polynomials v4.0.11
  [2913bbd2] StatsBase v0.34.3
  [d6d074c3] VideoIO v1.1.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [44cfe95a] Pkg v1.11.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.