diff --git a/Project.toml b/Project.toml index 92ec7a9..bb0dfb6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaveletsExt" uuid = "8f464e1e-25db-479f-b0a5-b7680379e03f" authors = ["Zeng Fung Liew ", "Shozen Dan "] -version = "0.1.14" +version = "0.1.18" [deps] AverageShiftedHistograms = "77b51b56-6f8f-5c3a-9cb4-d71f9594ea6e" @@ -13,6 +13,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" diff --git a/README.md b/README.md index b2311da..e7c4025 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,7 @@ [![deps](https://juliahub.com/docs/WaveletsExt/deps.svg)](https://juliahub.com/ui/Packages/WaveletsExt/iZ29j?t=2) [![version](https://juliahub.com/docs/WaveletsExt/version.svg)](https://juliahub.com/ui/Packages/WaveletsExt/iZ29j) [![pkgeval](https://juliahub.com/docs/WaveletsExt/pkgeval.svg)](https://juliahub.com/ui/Packages/WaveletsExt/iZ29j) +[![WaveletsExt Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/WaveletsExt)](https://pkgs.genieframework.com?packages=WaveletsExt) This package is a [Julia](https://github.com/JuliaLang/julia) extension package to [Wavelets.jl](https://github.com/JuliaDSP/Wavelets.jl) (WaveletsExt is short for Wavelets diff --git a/docs/Project.toml b/docs/Project.toml index 18276ce..929ea6a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,12 +1,16 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" WaveletsExt = "8f464e1e-25db-479f-b0a5-b7680379e03f" [compat] +BenchmarkTools = "1.2" Documenter = "0.27.6" Images = "0.24, 0.25" Plots = "1.21.3" diff --git a/docs/make.jl b/docs/make.jl index b7e4abd..ad850fd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,13 +13,15 @@ makedocs( "Transforms" => "manual/transforms.md", "Best Basis" => "manual/bestbasis.md", "Denoising" => "manual/denoising.md", - "Local Discriminant Basis" => "manual/localdiscriminantbasis.md" + "Local Discriminant Basis" => "manual/localdiscriminantbasis.md", + "Wavelets and Fast Numerical Algorithms" => "manual/wavemult.md" ], "API" => Any[ "DWT" => "api/dwt.md", "ACWT" => "api/acwt.md", "SWT" => "api/swt.md", "SIWPD" => "api/siwpd.md", + "WaveMult" => "api/wavemult.md", "Best Basis" => "api/bestbasis.md", "Denoising" => "api/denoising.md", "LDB" => "api/ldb.md", diff --git a/docs/src/api/wavemult.md b/docs/src/api/wavemult.md new file mode 100644 index 0000000..3152cf9 --- /dev/null +++ b/docs/src/api/wavemult.md @@ -0,0 +1,30 @@ +```@index +Modules = [WaveMult] +``` + +## Wavelet Multiplication +```@docs +WaveMult.nonstd_wavemult +WaveMult.std_wavemult +``` + +## Matrix to Sparse Format +```@docs +WaveMult.mat2sparseform_nonstd +WaveMult.mat2sparseform_std +``` + +## Fast Transform Algorithms +```@docs +WaveMult.ns_dwt +WaveMult.ns_idwt +WaveMult.sft +WaveMult.isft +``` + +## Miscellaneous +```@docs +WaveMult.dyadlength +WaveMult.ndyad +WaveMult.stretchmatrix +``` \ No newline at end of file diff --git a/docs/src/manual/wavemult.md b/docs/src/manual/wavemult.md new file mode 100644 index 0000000..0812b64 --- /dev/null +++ b/docs/src/manual/wavemult.md @@ -0,0 +1,84 @@ +# [Fast Numerical Algorithms using Wavelet Transforms](@id wavemult_manual) +This demonstration is based on the paper *[Fast wavelet transforms and numerical algorithms](https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160440202)* by G. Beylkin, R. Coifman, and V. Rokhlin and the lecture notes *[Wavelets and Fast Numerical ALgorithms](https://arxiv.org/abs/comp-gas/9304004)* by G. Beylkin. + +Assume we have a matrix ``M \in \mathbb{R}^{n \times n}``, and a vector ``x \in \mathbb{R}^n``, there are two ways to compute ``y = Mx``: + +1) Use the regular matrix multiplication, which takes ``O(n^2)`` time. +2) Transform both ``M`` and ``x`` to their standard/nonstandard forms ``\tilde M`` and ``\tilde x`` respectively. Compute the multiplication of ``\tilde y = \tilde M \tilde x``, then find the inverse of ``\tilde y``.If the matrix ``\tilde M`` is sparse, this can take ``O(n)`` time. + +## Examples and Comparisons +In our following examples, we compare the methods described in the previous section. We focus on two important metrics: time taken and accuracy of the algorithm. + +We start off by setting up the matrix `M`, the vector `x`, and the wavelet `wt`. As this algorithm is more efficient when ``n`` is large, we set ``n=2048``. +```@example wavemult +using Wavelets, WaveletsExt, Random, LinearAlgebra, BenchmarkTools + +# Construct matrix M, vector x, and wavelet wt +function data_setup(n::Integer, random_state = 1234) + M = zeros(n,n) + for i in axes(M,1) + for j in axes(M,2) + M[i,j] = i==j ? 0 : 1/abs(i-j) + end + end + + Random.seed!(random_state) + x = randn(n) + + return M, x +end + +M, x = data_setup(2048) +wt = wavelet(WT.haar) +nothing # hide +``` + +The following code block computes the regular matrix multiplication. +```@example wavemult +y₀ = M * x +nothing # hide +``` + +The following code block computes the nonstandard wavelet multiplication. Note that we will only be running 4 levels of wavelet transform, as running too many levels result in large computation time, without much improvement in accuracy. +```@example wavemult +L = 4 +NM = mat2sparseform_nonstd(M, wt, L) +y₁ = nonstd_wavemult(NM, x, wt, L) +nothing # hide +``` + +The following code block computes the standard wavelet multiplication. Once again, we will only be running 4 levels of wavelet transform here. +```@example wavemult +SM = mat2sparseform_std(M, wt, L) +y₂ = std_wavemult(SM, x, wt, L) +nothing # hide +``` + +!!! tip "Performance Tip" + Running more levels of transform (i.e. setting large values of `L`) does not necessarily result in more accurate approximations. Setting `L` to be a reasonably large value (e.g. 3 or 4) not only reduces computational time, but could also produce better results. + +We assume that the regular matrix multiplication produces the true results (bar some negligible machine errors), and we compare our results from the nonstandard and standard wavelet multiplications based on their relative errors (``\frac{||\hat y - y_0||_2}{||y_0||_2}``). +```@repl wavemult +relativenorm(y₁, y₀) # Comparing nonstandard wavelet multiplication with true value +relativenorm(y₂, y₀) # Comparing standard wavelet multiplication with true value +``` + +Lastly, we compare the computation time for each algorithm. +```@repl wavemult +@benchmark $M * $x # Benchmark regular matrix multiplication +@benchmark nonstd_wavemult($NM, $x, $wt, L) # Benchmark nonstandard wavelet multiplication +@benchmark std_wavemult($SM, $x, $wt, L) # Benchmark standard wavelet multiplication +``` + +As we can see above, the nonstandard and standard wavelet multiplication methods are significantly faster than the regular method and provides reasonable approximations to the true values. + +!!! tip "Performance Tips" + This method should only be used when the dimensions of ``M`` and ``x`` are large. Additionally, this would be more useful when one aims to compute ``Mx_0, Mx_1, \ldots, Mx_k`` where ``M`` remains constant throughout the ``k`` computations. In this case, it is more efficient to use + ```julia + NM = mat2sparseform_nonstd(M, wt) + y = nonstd_wavemult(NM, x, wt) + ``` + instead of + ```julia + y = nonstd_wavemult(M, x, wt) + ``` \ No newline at end of file diff --git a/paper/bestbasis.png b/paper/bestbasis.png index c415c32..064dfa9 100644 Binary files a/paper/bestbasis.png and b/paper/bestbasis.png differ diff --git a/paper/denoising.png b/paper/denoising.png index 7dc19c2..16e0c27 100644 Binary files a/paper/denoising.png and b/paper/denoising.png differ diff --git a/paper/examples.jl b/paper/examples.jl index e3b29eb..2720432 100644 --- a/paper/examples.jl +++ b/paper/examples.jl @@ -24,7 +24,7 @@ p2 = wiggle(y) |> p -> plot!(p, yticks=1:9, title="Stationary WT") # Combine and save plot p = plot(p1, p2, layout=(1,2)) -savefig(p, "transforms.png") +savefig(p, "tpaper/ransforms.png") # 2. Best Basis Algorithms ----------------------------------------------------------------- # Compared to Wavelets.jl, WaveletsExt.jl has an extended ability of catering toward @@ -39,12 +39,12 @@ xw = wpdall(x, wt, 6) # ----- Joint Best Basis (JBB) tree = bestbasistree(xw, JBB()) -p1 = plot_tfbdry(tree, 6, nd_col=:green, ln_col=:black, bg_col=:white) |> +p1 = plot_tfbdry(tree, 6, node_color=:green, line_color=:black, background_color=:white) |> p -> plot!(p, title="JBB") # ----- Least Statistically Dependent Basis (LSDB) tree = bestbasistree(xw, LSDB()) -p2 = plot_tfbdry(tree, 6, nd_col=:green, ln_col=:black, bg_col=:white) |> +p2 = plot_tfbdry(tree, 6, node_color=:green, line_color=:black, background_color=:white) |> p -> plot!(p, title="LSDB") # Combine and save plot @@ -118,7 +118,7 @@ ldb = LocalDiscriminantBasis(wt=wt, X̂ = fit_transform(ldb, X, y) # Plot the best basis for feature extraction -p2 = plot_tfbdry(ldb.tree, 6, nd_col=:green, ln_col=:black, bg_col=:white) +p2 = plot_tfbdry(ldb.tree, 6, node_color=:green, line_color=:black, background_color=:white) plot!(p2, title="Basis Selection using LDB") p = plot(p1, p2, size=(600,300)) diff --git a/paper/ldb.png b/paper/ldb.png index 70d18f8..baaa3e0 100644 Binary files a/paper/ldb.png and b/paper/ldb.png differ diff --git a/paper/paper.md b/paper/paper.md index 10e34cd..cd0a391 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -33,7 +33,7 @@ One of the most distinguishing features of `WaveletsExt.jl` is the presence of a # Examples ## 1. Redundant Wavelet Transforms -`WaveletsExt.jl` implements several redundant wavelet transforms including Stationary Wavelet Transform (SWT) [@Nason:1995] and Autocorrelation Wavelet Transform [@Saito:1993]. These transformations can be performed using the `acdwt` and `sdwt` functions, and the resulting decomposition can be visualized with the `wiggle` function included in `WaveletsExt.jl`. +`WaveletsExt.jl` implements several redundant wavelet transforms including Autocorrelation Wavelet Transform [@Saito:1993] and Stationary Wavelet Transform (SWT) [@Nason:1995]. These transformations can be performed using the `acdwt` and `sdwt` functions, and the resulting decomposition can be visualized with the `wiggle` function included in `WaveletsExt.jl`. ```julia using Plots, Wavelets, WaveletsExt @@ -56,7 +56,7 @@ savefig(p, "transforms.png") ``` !["Wiggle" plots displaying the value of coefficients at each level of the autocorrelation and stationary wavelet transform for a unit impulse signal. \label{fig:transforms}](transforms.png) -## Best Basis Algorithms +## 2. Best Basis Algorithms `WaveletsExt.jl` can select a best basis for a multiple signal input (i.e., an array of signals) through the Joint Best Basis (JBB) [@Wickerhauser:1996] or Least Statistically Dependent Basis (LSDB) [@Saito:2001] algorithms. The resulting best basis tree can be visualized using `plot_tfbdry` also included in `WaveletsExt.jl`. ```julia @@ -71,12 +71,18 @@ xw = wpdall(x, wt, 6) # ----- Joint Best Basis (JBB) tree = bestbasistree(xw, JBB()) -p1 = plot_tfbdry(tree, nd_col=:black, ln_col=:black, bg_col=:white) |> +p1 = plot_tfbdry(tree, + node_color=:green, + line_color=:black, + background_color=:white) |> p -> plot!(p, title="JBB") # ----- Least Statistically Dependent Basis (LSDB) tree = bestbasistree(xw, LSDB()) -p2 = plot_tfbdry(tree, nd_col=:black, ln_col=:black, bg_col=:white) |> +p2 = plot_tfbdry(tree, + node_color=:green, + line_color=:black, + background_color=:white) |> p -> plot!(p, title="LSDB") # Combine and save plot @@ -85,7 +91,7 @@ savefig(p, "bestbasis.png") ``` ![The best basis trees of 100 HeaviSine signals (A sinusoid + two Heaviside step functions) [@Donoho:1995a; @Donoho:1995b] selected by the JBB and LSDB algorithms. Each row represents a decomposition level, where level 0 is the original input signal, and each cell represents a frequency subband (low to high frequency from left to right). The colored cells indicate those subbands selected by the JBB (left) and the LSDB (right) algorithms. \label{fig:bestbasis}](bestbasis.png) -## Denoising Algorithms +## 3. Denoising Algorithms `WaveletsExt.jl` contains two functions for denoising: `denoise` and `denoiseall`. The former denoises a single signal input whereas the latter denoises multiple signal input. For more examples of denoising algorithms in `WaveletsExt.jl`, see [@Liew:2021]. ```julia @@ -126,7 +132,7 @@ savefig(p, "denoising.png") ``` ![Left: HeaviSine signals with Gaussian noise. Black lines represent the original (non-noisy) signal. Right: Simultaneously denoised signals using the JBB algorithm with a universal thresholding constant determined by the VisuShrink method [@Donoho:1994]. \label{fig:denoising}](denoising.png) -## Feature Extraction +## 4. Feature Extraction For signal classification problems, users can extract distinguishing features localized in both the time and frequency domains using the Local Discriminant Basis (LDB) algorithm. Further details can be found in the original papers by Saito and his collaborators [@Saito:1995; @Saito:2002] as well as the interactive tutorial [@Dan:2021]. ```julia @@ -159,7 +165,10 @@ ldb = LocalDiscriminantBasis( X̂ = fit_transform(ldb, X, y) # Plot the best basis for feature extraction -p2 = plot_tfbdry(ldb.tree, nd_col=:black, ln_col=:black, bg_col=:white) +p2 = plot_tfbdry(ldb.tree, + node_color=:green, + line_color=:black, + background_color=:white) plot!(p2, title="Basis Selection using LDB") # Combine and save plot diff --git a/src/WaveletsExt.jl b/src/WaveletsExt.jl index 38d7f4b..429e5dc 100644 --- a/src/WaveletsExt.jl +++ b/src/WaveletsExt.jl @@ -11,6 +11,7 @@ include("mod/SWT.jl") include("mod/Denoising.jl") include("mod/LDB.jl") include("mod/Visualizations.jl") +include("mod/WaveMult.jl") using Reexport @reexport using .DWT, @@ -21,6 +22,7 @@ using Reexport .SWT, .SIWPD, .ACWT, - .Visualizations + .Visualizations, + .WaveMult end \ No newline at end of file diff --git a/src/mod/ACWT.jl b/src/mod/ACWT.jl index 9c9ff71..a14d196 100644 --- a/src/mod/ACWT.jl +++ b/src/mod/ACWT.jl @@ -1006,8 +1006,8 @@ end # return (v₀ + v₁) / √2 # end -include("acwt_utils.jl") -include("acwt_one_level.jl") -include("acwt_all.jl") +include("acwt/acwt_utils.jl") +include("acwt/acwt_one_level.jl") +include("acwt/acwt_all.jl") end # end module \ No newline at end of file diff --git a/src/mod/BestBasis.jl b/src/mod/BestBasis.jl index efa9af7..d0fc3f2 100644 --- a/src/mod/BestBasis.jl +++ b/src/mod/BestBasis.jl @@ -33,8 +33,8 @@ using ..DWT, ..SIWPD -include("bestbasis_costs.jl") -include("bestbasis_tree.jl") +include("bestbasis/bestbasis_costs.jl") +include("bestbasis/bestbasis_tree.jl") ## ----- BEST TREE SELECTION ----- # Tree selection for 1D signals diff --git a/src/mod/DWT.jl b/src/mod/DWT.jl index 246016c..97a1abf 100644 --- a/src/mod/DWT.jl +++ b/src/mod/DWT.jl @@ -709,7 +709,7 @@ function Wavelets.Transforms.iwpt!(x̂::AbstractArray{T,2}, return x̂ end -include("dwt_one_level.jl") -include("dwt_all.jl") +include("dwt/dwt_one_level.jl") +include("dwt/dwt_all.jl") end # end module \ No newline at end of file diff --git a/src/mod/LDB.jl b/src/mod/LDB.jl index ad5825a..1e6c308 100644 --- a/src/mod/LDB.jl +++ b/src/mod/LDB.jl @@ -41,8 +41,8 @@ using ..Utils, ..DWT import ..BestBasis: bestbasis_treeselection -include("ldb_energymap.jl") -include("ldb_measures.jl") +include("ldb/ldb_energymap.jl") +include("ldb/ldb_measures.jl") ## LOCAL DISCRIMINANT BASIS """ diff --git a/src/mod/SWT.jl b/src/mod/SWT.jl index 2b39fad..ca0e67d 100644 --- a/src/mod/SWT.jl +++ b/src/mod/SWT.jl @@ -1198,7 +1198,7 @@ function iswpd!(x::AbstractMatrix{T}, return x end -include("swt_one_level.jl") -include("swt_all.jl") +include("swt/swt_one_level.jl") +include("swt/swt_all.jl") end # end module \ No newline at end of file diff --git a/src/mod/Utils.jl b/src/mod/Utils.jl index 52f9161..5f319a4 100644 --- a/src/mod/Utils.jl +++ b/src/mod/Utils.jl @@ -541,8 +541,8 @@ function getcolrange(n::Integer, idx::T) where T<:Integer end end -include("utils_tree.jl") -include("utils_metrics.jl") -include("utils_dataset.jl") +include("utils/utils_tree.jl") +include("utils/utils_metrics.jl") +include("utils/utils_dataset.jl") end # end module \ No newline at end of file diff --git a/src/mod/WaveMult.jl b/src/mod/WaveMult.jl new file mode 100644 index 0000000..6f21192 --- /dev/null +++ b/src/mod/WaveMult.jl @@ -0,0 +1,19 @@ +module WaveMult +export mat2sparseform_std, + mat2sparseform_nonstd, + ns_dwt, + ns_idwt, + std_wavemult, + nonstd_wavemult + +using Wavelets, + LinearAlgebra, + SparseArrays + +import ..DWT: dwt_step!, idwt_step! + +include("wavemult/utils.jl") +include("wavemult/mat2sparse.jl") +include("wavemult/wavemult.jl") +include("wavemult/transforms.jl") +end # module \ No newline at end of file diff --git a/src/mod/acwt_all.jl b/src/mod/acwt/acwt_all.jl similarity index 100% rename from src/mod/acwt_all.jl rename to src/mod/acwt/acwt_all.jl diff --git a/src/mod/acwt_one_level.jl b/src/mod/acwt/acwt_one_level.jl similarity index 100% rename from src/mod/acwt_one_level.jl rename to src/mod/acwt/acwt_one_level.jl diff --git a/src/mod/acwt_utils.jl b/src/mod/acwt/acwt_utils.jl similarity index 100% rename from src/mod/acwt_utils.jl rename to src/mod/acwt/acwt_utils.jl diff --git a/src/mod/bestbasis_costs.jl b/src/mod/bestbasis/bestbasis_costs.jl similarity index 100% rename from src/mod/bestbasis_costs.jl rename to src/mod/bestbasis/bestbasis_costs.jl diff --git a/src/mod/bestbasis_tree.jl b/src/mod/bestbasis/bestbasis_tree.jl similarity index 100% rename from src/mod/bestbasis_tree.jl rename to src/mod/bestbasis/bestbasis_tree.jl diff --git a/src/mod/dwt_all.jl b/src/mod/dwt/dwt_all.jl similarity index 100% rename from src/mod/dwt_all.jl rename to src/mod/dwt/dwt_all.jl diff --git a/src/mod/dwt_one_level.jl b/src/mod/dwt/dwt_one_level.jl similarity index 100% rename from src/mod/dwt_one_level.jl rename to src/mod/dwt/dwt_one_level.jl diff --git a/src/mod/ldb_energymap.jl b/src/mod/ldb/ldb_energymap.jl similarity index 100% rename from src/mod/ldb_energymap.jl rename to src/mod/ldb/ldb_energymap.jl diff --git a/src/mod/ldb_measures.jl b/src/mod/ldb/ldb_measures.jl similarity index 100% rename from src/mod/ldb_measures.jl rename to src/mod/ldb/ldb_measures.jl diff --git a/src/mod/swt_all.jl b/src/mod/swt/swt_all.jl similarity index 100% rename from src/mod/swt_all.jl rename to src/mod/swt/swt_all.jl diff --git a/src/mod/swt_one_level.jl b/src/mod/swt/swt_one_level.jl similarity index 100% rename from src/mod/swt_one_level.jl rename to src/mod/swt/swt_one_level.jl diff --git a/src/mod/utils_dataset.jl b/src/mod/utils/utils_dataset.jl similarity index 100% rename from src/mod/utils_dataset.jl rename to src/mod/utils/utils_dataset.jl diff --git a/src/mod/utils_metrics.jl b/src/mod/utils/utils_metrics.jl similarity index 100% rename from src/mod/utils_metrics.jl rename to src/mod/utils/utils_metrics.jl diff --git a/src/mod/utils_tree.jl b/src/mod/utils/utils_tree.jl similarity index 100% rename from src/mod/utils_tree.jl rename to src/mod/utils/utils_tree.jl diff --git a/src/mod/wavemult/mat2sparse.jl b/src/mod/wavemult/mat2sparse.jl new file mode 100644 index 0000000..f54587f --- /dev/null +++ b/src/mod/wavemult/mat2sparse.jl @@ -0,0 +1,100 @@ +""" + mat2sparseform_nonstd(M, wt, [L], [ϵ]) + +Transform the matrix `M` into the wavelet basis. Then, it is stretched into its nonstandard +form. Elements exceeding `ϵ * maximum column norm` are set to zero. The resulting output +sparse matrix, and can be used as input to `nonstd_wavemult`. + +# Arguments +- `M::AbstractMatrix{T} where T<:AbstractFloat`: `n` by `n` matrix (`n` dyadic) to be put in + Sparse Nonstandard form. +- `wt::OrthoFilter`: Wavelet filter. +- `L::Integer`: (Default: `maxtransformlevels(M)`) Number of decomposition levels. +- `ϵ::T where T<:AbstractFloat`: (Default: `1e-4`) Truncation Criterion. + +# Returns +- `NM::SparseMatrixCSC{T, Integer}`: Sparse nonstandard form of matrix of size 2n x 2n. + +# Examples +```jldoctest +julia> using Wavelets, WaveletsExt; import Random: seed! + +julia> seed!(1234); M = randn(4,4); wt = wavelet(WT.haar); + +julia> mat2sparseform_nonstd(M, wt) +8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 16 stored entries: + 1.88685 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ 0.363656 ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ 2.49634 -1.08139 ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0187 0.539411 + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.68141 0.0351839 + ⋅ ⋅ ⋅ ⋅ -1.39713 -1.21352 0.552745 0.427717 + ⋅ ⋅ ⋅ ⋅ -1.05882 0.16666 -0.124156 -0.218902 +``` + +**See also:** [`mat2sparseform_std`](@ref), [`stretchmatrix`](@ref) +""" +function mat2sparseform_nonstd(M::AbstractMatrix{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(M), + ϵ::T = 1e-4) where T<:AbstractFloat + @assert size(M,1) == size(M,2) + n = size(M, 1) + Mw = dwt(M, wt, L) # 2D wavelet transform + # Find maximum column norm of Mw + maxcolnorm = mapslices(norm, Mw, dims = 1) |> maximum + nilMw = Mw .* (abs.(Mw) .> ϵ * maxcolnorm) # "Remove" values close to zero + nz_vals = nilMw[nilMw .≠ 0] # Extract all nonzero values + nz_indx = findall(!iszero, nilMw) # Extract indices of nonzero values + i = getindex.(nz_indx, 1) + j = getindex.(nz_indx, 2) + ie, je = stretchmatrix(i, j, n, L) # Stretch matrix Mw + NM = sparse(ie, je, nz_vals, 2*n, 2*n) # Convert to sparse matrix + return NM +end + +""" + mat2sparseform_std(M, wt, [L], [ϵ]) + +Transform the matrix `M` into the standard form. Then, elements exceeding `ϵ * maximum +column norm` are set to zero. The resulting output sparse matrix, and can be used as input +to `std_wavemult`. + +# Arguments +- `M::AbstractMatrix{T} where T<:AbstractFloat`: Matrix to be put in Sparse Standard form. +- `wt::OrthoFilter`: Wavelet filter. +- `L::Integer`: (Default: `maxtransformlevels(M)`) Number of decomposition levels. +- `ϵ::T where T<:AbstractFloat`: (Default: `1e-4`) Truncation Criterion. + +# Returns +- `SM::SparseMatrixCSC{T, Integer}`: Sparse standard form of matrix of size ``n \times n``. + +# Examples +```jldoctest +julia> using Wavelets, WaveletsExt; import Random: seed! + +julia> seed!(1234); M = randn(4,4); wt = wavelet(WT.haar); + +julia> mat2sparseform_std(M, wt) +4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 16 stored entries: + 1.88685 0.363656 0.468602 0.4063 + 2.49634 -1.08139 1.90927 -0.356542 + -1.84601 0.129829 0.552745 0.427717 + -0.630852 0.866545 -0.124156 -0.218902 +``` + +**See also:** [`mat2sparseform_nonstd`](@ref), [`sft`](@ref) +""" +function mat2sparseform_std(M::AbstractMatrix{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(M), + ϵ::T = 1e-4) where T<:AbstractFloat + @assert size(M,1) == size(M,2) + Mw = sft(M, wt, L) # Transform to standard form + # Find maximum column norm of Mw + maxcolnorm = mapslices(norm, Mw, dims = 1) |> maximum + nilMw = Mw .* (abs.(Mw) .> ϵ * maxcolnorm) # Remove values close to zero + SM = sparse(nilMw) # Convert to sparse matrix + return SM +end \ No newline at end of file diff --git a/src/mod/wavemult/transforms.jl b/src/mod/wavemult/transforms.jl new file mode 100644 index 0000000..e7e16d5 --- /dev/null +++ b/src/mod/wavemult/transforms.jl @@ -0,0 +1,228 @@ +""" + ns_dwt(x, wt, [L]) + +Nonstandard wavelet transform on 1D signals. + +# Arguments +- `x::AbstractVector{T} where T<:AbstractFloat`: 1D signal of length ``2^J``. +- `wt::OrthoFilter`: Type of wavelet filter. +- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of decomposition levels. + +# Returns +- `nxw::Vector{T}`: Nonstandard wavelet transform of `x` of length ``2^{J+1}``. + +# Examples +```jldoctest +julia> using Wavelets, WaveletsExt; import Random: seed! + +julia> seed!(1234); x = randn(4) +4-element Vector{Float64}: + -0.3597289068234817 + 1.0872084924285859 + -0.4195896169388487 + 0.7189099374659392 + +julia> wt = wavelet(WT.haar); + +julia> nxw = ns_dwt(x, wt) # Nonstandard transform +8-element Vector{Float64}: + 0.0 + 0.0 + 0.5133999530660974 + -0.2140796325390069 + 0.5144057481561487 + 0.21165142839163664 + 1.0231392469635638 + 0.8050407552974883 + +julia> x̂ = ns_idwt(nxw, wt) # Nonstandard inverse transform +4-element Vector{Float64}: + 0.004010885979070622 + 1.4509482852311382 + -0.2699294566753035 + 0.8685700977294846 + +julia> x ≈ x̂ # Unlike standard dwt, x ≠ x̂ +false +``` + +**See also:** [`nonstd_wavemult`](@ref), [`mat2sparseform_nonstd`](@ref), [`ns_idwt`](@ref), +[`ndyad`](@ref) +""" +function ns_dwt(x::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(x)) where T<:AbstractFloat + Lmax = maxtransformlevels(x) + n = length(x) + @assert 1 ≤ L ≤ Lmax + @assert ispow2(n) + + nxw = zeros(T, 2*n) + g, h = WT.makereverseqmfpair(wt, true) + for l in 1:L + w₁ = @view nxw[ndyad(l, Lmax, false)] + w₂ = @view nxw[ndyad(l, Lmax, true)] + v = l == 1 ? x : @view nxw[ndyad(l-1, Lmax, false)] + dwt_step!(w₁, w₂, v, h, g) + end + nxw[1:1<<(Lmax-L)] = nxw[ndyad(L, Lmax, false)] + return nxw +end + +""" + ns_idwt(nxw, wt, [L]) + +Inverse nonstandard wavelet transform on 1D signals. + +# Arguments +- `nxw::AbstractVector{T} where T<:AbstractFloat`: Nonstandard wavelet transformed 1D signal + of length ``2^{J+1}``. +- `wt::OrthoFilter`: Type of wavelet filter. +- `L::Integer`: (Default: `maxtransformlevels(nxw) - 1`) Number of decomposition levels. + +# Returns +- `x::Vector{T}`: 1D signal of length ``2^J``. + +# Examples +```jldoctest +julia> using Wavelets, WaveletsExt; import Random: seed! + +julia> seed!(1234); x = randn(4) +4-element Vector{Float64}: + -0.3597289068234817 + 1.0872084924285859 + -0.4195896169388487 + 0.7189099374659392 + +julia> wt = wavelet(WT.haar); + +julia> nxw = ns_dwt(x, wt) # Nonstandard transform +8-element Vector{Float64}: + 0.0 + 0.0 + 0.5133999530660974 + -0.2140796325390069 + 0.5144057481561487 + 0.21165142839163664 + 1.0231392469635638 + 0.8050407552974883 + +julia> x̂ = ns_idwt(nxw, wt) # Nonstandard inverse transform +4-element Vector{Float64}: + 0.004010885979070622 + 1.4509482852311382 + -0.2699294566753035 + 0.8685700977294846 + +julia> x ≈ x̂ # Unlike standard dwt, x ≠ x̂ +false +``` + +**See also:** [`nonstd_wavemult`](@ref), [`mat2sparseform_nonstd`](@ref), [`ns_dwt`](@ref), +[`ndyad`](@ref) +""" +function ns_idwt(nxw::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(nxw)-1) where T<:AbstractFloat + Lmax = maxtransformlevels(nxw) - 1 + n = length(nxw) ÷ 2 + @assert 1 ≤ L ≤ Lmax + @assert ispow2(n) + + x = zeros(T, n) + x[1:1<<(Lmax-L)] = nxw[1:1<<(Lmax-L)] + g, h = WT.makereverseqmfpair(wt, true) + for l in L:-1:1 + w₁ = nxw[ndyad(l, Lmax, false)] + x[1:1<<(Lmax-l)] + w₂ = @view nxw[ndyad(l, Lmax, true)] + v = @view x[1:1<<(Lmax-l+1)] + idwt_step!(v, w₁, w₂, h, g) + end + return x +end + +""" + sft(M, wt, [L]) + +Transforms a matrix `M` to be then represented in the sparse standard form. This is achieved +by first computing ``L`` levels of wavelet transform on each column of `M`, and then +computing ``L`` levels of wavelet transform on each row of `M`. + +# Arguments +- `M::AbstractMatrix{T} where T<:AbstractFloat`: Matrix to be put in standard form. +- `wt::OrthoFilter`: Type of wavelet filter. +- `L::Integer`: Number of decomposition levels. + +# Returns +- `Mw::Matrix{T}`: Matrix in the standard form. + +# Examples +```julia +import WaveletsExt.WaveMult: sft +using Wavelets + +M = randn(4,4); wt = wavelet(WT.haar) +Mw = sft(M, wt) +M̂ = isft(Mw, wt) +``` + +**See also:** [`mat2sparseform_std`](@ref), [`isft`](@ref) +""" +function sft(M::AbstractMatrix{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(M)) where T<:AbstractFloat + @assert 1 ≤ L ≤ maxtransformlevels(M) + n, m = size(M) + + Mw = similar(M) + for j in 1:m # Transform each column + Mw[:,j] = dwt(M[:,j], wt, L) + end + for i in 1:n # Transform each row + Mw[i,:] = dwt(Mw[i,:], wt, L) + end + return Mw +end + +""" + isft(Mw, wt, [L]) + +Reconstructs the matrix `M` from the sparse standard form `Mw`. This is achieved by first +computing ``L`` levels of inverse wavelet transform on each row of `Mw`, and then computing +``L`` levels of inverse wavelet transform on each column of `Mw`. + +# Arguments +- `Mw::AbstractMatrix{T} where T<:AbstractFloat`: Matrix in standard form. +- `wt::OrthoFilter`: Type of wavelet filter. +- `L::Integer`: Number of decomposition levels. + +# Returns +- `M::Matrix{T}`: Reconstructed matrix. + +# Examples +```julia +import WaveletsExt.WaveMult: sft +using Wavelets + +M = randn(4,4); wt = wavelet(WT.haar) +Mw = sft(M, wt) +M̂ = isft(Mw, wt) +``` + +**See also:** [`sft`](@ref) +""" +function isft(Mw::AbstractMatrix{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(Mw)) where T<:AbstractFloat + @assert 1 ≤ L ≤ maxtransformlevels(Mw) + n, m = size(Mw) + + M = similar(Mw) + for i in 1:n # Transform each row + M[i,:] = idwt(Mw[i,:], wt, L) + end + for j in 1:m # Transform each column + M[:,j] = idwt(M[:,j], wt, L) + end + return M +end \ No newline at end of file diff --git a/src/mod/wavemult/utils.jl b/src/mod/wavemult/utils.jl new file mode 100644 index 0000000..5ebb450 --- /dev/null +++ b/src/mod/wavemult/utils.jl @@ -0,0 +1,155 @@ +""" + dyadlength(x) + dyadlength(n) + +Find dyadic length of array. + +# Arguments +- `x::AbstractArray`: Array of length `n`. Preferred array length is ``2^J`` where ``J`` is + an integer. +- `n::Integer`: Length of array. + +# Returns +- `J::Integer`: Least power of two greater than `n`. + +!!! note + The function `dyadlength` is very similar to the function `maxtransformlevels` from + Wavelets.jl. The only difference here is the way it handles `n` when `n` is not a power + of 2. The example below provides a demonstration of the differences in the 2 functions. + +# Examples +```jldoctest +julia> import WaveletsExt.WaveMult: dyadlength + +julia> import Wavelets: maxtransformlevels + +julia> x = randn(16); dyadlength(x) +4 + +julia> dyadlength(16) # Same as previous +4 + +julia> maxtransformlevels(16) # Equivalent to dyadlength when n is power of 2 +4 + +julia> dyadlength(15) # Produces warning when n is not power of 2 +┌ Warning: Dyadlength n != 2^J +└ @ WaveletsExt.WaveMult +4 + +julia> maxtransformlevels(15) # Not equivalent to dyadlength when n is not power of 2 +0 +``` +""" +function dyadlength(n::T) where T<:Integer + J = ceil(T, log2(n)) + if 1< M = [1 0 0 0; + 0 2 0 0; + 0 0 3 0; + 0 0 0 4]; + +julia> idx = findall(!iszero, M) +2-element Vector{CartesianIndex{2}}: + CartesianIndex(1, 1) + CartesianIndex(2, 2) + +julia> i = getindex.(idx, 1) +4-element Vector{Int64}: + 1 + 2 + 3 + 4 + +julia> j = getindex.(idx, 2) +4-element Vector{Int64}: + 1 + 2 + 3 + 4 + +julia> stretchmatrix(i, j, 4, 2) +([1, 4, 7, 8], [1, 4, 7, 8]) +``` +""" +function stretchmatrix(i::AbstractVector{T}, j::AbstractVector{T}, n::T, L::T) where + T<:Integer + Lmax = maxtransformlevels(n) + @assert 1 ≤ L ≤ Lmax + ie = copy(i) + je = copy(j) + for l in 0:(L-1) + k = Lmax - l - 1 + cond = ((ie .> 1< 1< import WaveletsExt.WaveMult: ndyad + +julia> ndyad(1, 4, false) +17:24 + +julia> ndyad(1, 4, true) +25:32 +``` +""" +function ndyad(L::T, Lmax::T, gender::Bool) where T<:Integer + @assert L ≤ Lmax # Current level cannot be larger than max level + @assert L ≥ 1 # Current level must be at least 1 + k = Lmax - L + if gender + return (1<<(k+1) + 1< using Wavelets, WaveletsExt + +julia> M = randn(4,4); x = randn(4); wt = wavelet(WT.haar); + +julia> NM = mat2sparseform_nonstd(M, wt); y₀ = nonstd_wavemult(NM, x, wt) +4-element Vector{Float64}: + -1.2590015844047044 + -1.3234024176418535 + 2.1158027198405627 + 1.5364835417087566 + +julia> y₁ = nonstd_wavemult(M, x, wt) +4-element Vector{Float64}: + -1.2590015844047044 + -1.3234024176418535 + 2.1158027198405627 + 1.5364835417087566 + +julia> y₀ == y₁ # Both methods are equivalent +true +``` + +**See also:* [`std_wavemult`](@ref), [`mat2sparseform_nonstd`](@ref), [`ns_dwt`](@ref), +[`ns_idwt`](@ref) +""" +function nonstd_wavemult(M::AbstractMatrix{T}, + x::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(x), + ϵ::T = 1e-4) where T<:AbstractFloat + NM = mat2sparseform_nonstd(M, wt, L, ϵ) + return nonstd_wavemult(NM, x, wt, L) +end + +function nonstd_wavemult(NM::SparseMatrixCSC{T,S}, + x::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(x)) where + {T<:AbstractFloat, S<:Integer} + nx = ns_dwt(x, wt, L) + ny = NM * nx + y = ns_idwt(ny, wt, L) + return y +end + +""" + std_wavemult(M, x, wt, [L], [ϵ]) + std_wavemult(SM, x, wt, [L]) + +If `M` is an ``n`` by ``n`` matrix, there are two ways to compute ``y = Mx``. The first is +to use the standard matrix multiplication to compute the product. This algorithm works in +order of ``O(n^2)``, where ``n`` is the length of `x`. + +The second is to transform both the matrix and vector to their standard forms and multiply +the standard forms. If the matrix is sparse in standard form, this can be an order ``O(n)`` +algorithm. + +!!! tip + One may choose to use the original matrix `M` as input by doing `std_wavemult(M, x, + wt, [L], [ϵ])`. However, if the standard form sparse matrix `SM` is already computed + prior, one can skip the redundant step by doing `std_wavemult(SM, x, wt, [L])`. + +# Arguments +- `M::AbstractVector{T} where T<:AbstractFloat`: ``n`` by ``n`` matrix. +- `NM::SparseMatrixCSC{T,S} where {T<:AbstractFloat, S<:Integer}`: Standard transformed + sparse matrix of `M`. +- `x::AbstractVector{T} where T<:AbstractFloat`: Vector of length ``n`` in natural basis. +- `wt::OrthoFilter`: Type of wavelet filter. +- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of decomposition levels. +- `ϵ::T where T<:AbstractFloat`: (Default: `1e-4`) Truncation criterion for standard + transform of `M`. + +# Returns +- `y::Vector{T}`: Standard form approximation of ``Mx``. + +# Examples +```jldoctest +julia> using Wavelets, WaveletsExt + +julia> M = randn(4,4); x = randn(4); wt = wavelet(WT.haar); + +julia> SM = mat2sparseform_std(M, wt); y₀ = std_wavemult(SM, x, wt) +4-element Vector{Float64}: + 2.2303830532617344 + -0.12704611958926648 + 2.656411941014368 + -4.811388406857621 + +julia> y₁ = std_wavemult(M, x, wt) +4-element Vector{Float64}: + 2.2303830532617344 + -0.12704611958926648 + 2.656411941014368 + -4.811388406857621 + +julia> y₀ == y₁ +true +``` + +**See also:** [`nonstd_wavemult`](@ref), [`mat2sparseform_std`](@ref) +""" +function std_wavemult(M::AbstractMatrix{T}, + x::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(x), + ϵ::T = 1e-4) where T<:AbstractFloat + SM = mat2sparseform_std(M, wt, L, ϵ) + return std_wavemult(SM, x, wt, L) +end + +function std_wavemult(SM::SparseMatrixCSC{T,S}, + x::AbstractVector{T}, + wt::OrthoFilter, + L::Integer = maxtransformlevels(x)) where + {T<:AbstractFloat, S<:Integer} + nx = dwt(x, wt, L) + ny = SM * nx + y = idwt(ny, wt, L) + return y +end \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index c65d7d1..6a416b7 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ImageQualityIndexes = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" diff --git a/test/runtests.jl b/test/runtests.jl index 2ab16d6..8395460 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,11 +9,13 @@ using Plots, Statistics, Wavelets, - WaveletsExt + WaveletsExt, + SparseArrays @testset "Utils" begin include("utils.jl") end -@testset "Visualizations" begin include("visualizations.jl") end @testset "Transforms" begin include("transforms.jl") end +@testset "Wavelet Multiplication" begin include("wavemult.jl") end @testset "Best Basis" begin include("bestbasis.jl") end @testset "Denoising" begin include("denoising.jl") end -@testset "LDB" begin include("ldb.jl") end \ No newline at end of file +@testset "LDB" begin include("ldb.jl") end +@testset "Visualizations" begin include("visualizations.jl") end \ No newline at end of file diff --git a/test/wavemult.jl b/test/wavemult.jl new file mode 100644 index 0000000..70a1b23 --- /dev/null +++ b/test/wavemult.jl @@ -0,0 +1,89 @@ +import WaveletsExt.WaveMult: dyadlength, stretchmatrix, ndyad, sft, isft + +@testset "Utilities" begin + # dyadlength + @test dyadlength(zeros(16)) == dyadlength(16) == 4 + @test_logs (:warn, "Dyadlength n != 2^J") dyadlength(15) + + # stretchmatrix + i = [1,2,3,4]; j = copy(i) + @test stretchmatrix(i,j,4,2) == ([1, 4, 7, 8], [1, 4, 7, 8]) + @test_throws AssertionError stretchmatrix(i,j,4,3) + @test_throws AssertionError stretchmatrix(i,j,4,0) + @test_throws AssertionError stretchmatrix(i,j,4,-1) + + # ndyad + @test ndyad(1,4,false) == 17:24 + @test ndyad(1,4,true) == 25:32 + @test_throws AssertionError ndyad(5,4,true) + @test_throws AssertionError ndyad(5,4,false) + @test_throws AssertionError ndyad(0,4,false) + @test_throws AssertionError ndyad(0,4,false) +end + +@testset "Wavelet Transforms" begin + # 1D Nonstandard transform + x = [1,2,-3,4.0] + wt = wavelet(WT.haar) + y = [2, 0, 2, -1, 2.1213, 0.7071, 0.7071, 4.9497] + z = [3.5, 4.5, -1.5, 5.5] + @test round.(ns_dwt(x, wt), digits = 4) == y + @test round.(ns_idwt(y, wt), digits = 4) == z + @test_throws AssertionError ns_dwt(x, wt, 3) + @test_throws AssertionError ns_dwt(x, wt, 0) + @test_throws AssertionError ns_idwt(y, wt, 3) + @test_throws AssertionError ns_idwt(y, wt, 0) + + # Standard form transform + x = [1 2 -3 4; 2 3 -4 1; 3 4 -1 2; 3 1 -2 3.0] + y = [4.75 -4.75 0.3536 7.0711; + 1.75 0.25 -1.0607 -1.4142; + -0.7071 -2.1213 0 -1; + -1.0607 1.0607 -1.5 1] + @test round.(sft(x, wt), digits = 4) == y + @test round.(isft(y, wt), digits = 4) == x + @test_throws AssertionError sft(x, wt, 3) + @test_throws AssertionError sft(x, wt, 0) + @test_throws AssertionError isft(y, wt, 3) + @test_throws AssertionError isft(y, wt, 0) + + # Build sparse matrices + x = [1 0 -3 0; 0 3 0 1; 3 0 -1 0; 0 1 0 3.0] + y = [2 0 0 0 0 0 0 0; + 0 0 0 0 0 0 0 0; + 0 0 0 -2 0 0 0 0; + 0 0 1 1 0 0 0 0; + 0 0 0 0 0 0 1 2; + 0 0 0 0 0 0 -1 2; + 0 0 0 0 1 2 2 -1; + 0 0 0 0 -1 2 2 1.0] + z = [2 -2 0 2.8284; + 1 1 -1.4142 0; + 2.1213 0.7071 2 -1; + 0.7071 2.1213 2 1] + y_sparse = sparse(y) + z_sparse = sparse(z) + @test round.(mat2sparseform_nonstd(x, wt), digits = 4) == y_sparse + @test round.(mat2sparseform_std(x, wt), digits = 4) == z_sparse + @test_throws AssertionError mat2sparseform_nonstd(randn(4,5), wt) + @test_throws AssertionError mat2sparseform_std(randn(4,5), wt) +end + +@testset "Wavelet Multiplication" begin + M = zeros(4,4) + for i in axes(M,1) + for j in axes(M,2) + M[i,j] = i==j ? 0 : 1/abs(i-j) + end + end + x = [0.5, 1, 1.5, 2] + wt = wavelet(WT.haar) + + # Test on Nonstandard Wavelet Multiplication + y = [2.4167, 3, 3.25, 2.1667] + @test round.(nonstd_wavemult(M, x, wt), digits = 4) == y + + # Test on Standard Wavelet Multiplication + y = [2.4167, 3, 3.25, 2.1667] + @test round.(std_wavemult(M, x, wt), digits = 4) == y +end \ No newline at end of file