Eye([domain_type=Float64::Type,] dim_in::Tuple)
Eye([domain_type=Float64::Type,] dims...)

Create the identity operator.

julia> op = Eye(Float64,(4,))
I  ℝ^4 -> ℝ^4

julia> op = Eye(2,3,4)
I  ℝ^(2, 3, 4) -> ℝ^(2, 3, 4)

julia> op*ones(2,3,4) == ones(2,3,4)
true
	

Zeros(domain_type::Type, dim_in::Tuple, [codomain_type::Type,] dim_out::Tuple)

Create a LinearOperator which, when multiplied with an array x of size dim_in, returns an array y of size dim_out filled with zeros.

For convenience Zeros can be constructed from any AbstractOperator.

julia> Zeros(Eye(10,20))
0  ℝ^(10, 20) -> ℝ^(10, 20)

julia> Zeros([Eye(10,20) Eye(10,20)])
[0,0]  ℝ^(10, 20)  ℝ^(10, 20) -> ℝ^(10, 20)
	

DiagOp(domain_type::Type, dim_in::Tuple, d::AbstractArray)
DiagOp(d::AbstractArray)

Creates a LinearOperator which, when multiplied with an array x, returns the elementwise product d.*x.

julia> D = DiagOp(Float64, (2, 2,), [1. 2.; 3. 4.])
╲  ℝ^(2, 2) -> ℝ^(2, 2)

julia> D*ones(2,2)
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0
	

GetIndex([domaintype=Float64::Type,] dimin::Tuple, idx...) GetIndex(x::AbstractArray, idx::Tuple)

Creates a LinearOperator which, when multiplied with x, returns x[idx]. Supported indices are:

• Colon(): selects the entire dimension.
• BitArray or AbstractArray{Bool}: selects elements where the mask is true.
• AbstractVector{Int}: selects elements by their indices.
• AbstractVector{CartesianIndex}: selects elements by their Cartesian indices.
• Tuple: a tuple of indices, where each element can be one of the above types.

julia> x = (1:10) .* 1.0;

julia> G = GetIndex(Float64,(10,), 1:3)
↓  ℝ^10 -> ℝ^3

julia> G*x
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> GetIndex(randn(10,20,30),(1:2,1:4,1))
↓  ℝ^(10, 20, 30) -> ℝ^(2, 4)

MatrixOp(domain_type=Float64::Type, dim_in::Tuple, A::AbstractMatrix)
MatrixOp(A::AbstractMatrix)
MatrixOp(A::AbstractMatrix, n_colons)

Creates a LinearOperator which, when multiplied with a vector x::AbstractVector, returns the product A*x.

The input x can be also a matrix: the number of columns must be given either in the second entry of dim_in::Tuple or using the constructor MatrixOp(A::AbstractMatrix, n_colons).

julia> MatrixOp(Float64,(10,),randn(20,10))
▒  ℝ^10 -> ℝ^20

julia> MatrixOp(randn(20,10))
▒  ℝ^10 -> ℝ^20

julia> MatrixOp(Float64,(10,20),randn(20,10))
▒  ℝ^(10, 20) -> ℝ^(20, 20)

julia> MatrixOp(randn(20,10),4)
▒  ℝ^(10, 4) -> ℝ^(20, 4)
	

LMatrixOp(domain_type=Float64::Type, dim_in::Tuple, b::Union{AbstractVector,AbstractMatrix})
LMatrixOp(b::AbstractVector, number_of_rows::Int)

Creates a LinearOperator which, when multiplied with a matrix X::AbstractMatrix, returns the product X*b.

julia> op = LMatrixOp(Float64,(3,4),ones(4))
(⋅)b  ℝ^(3, 4) -> ℝ^3

julia> op = LMatrixOp(ones(4),3)
(⋅)b  ℝ^(3, 4) -> ℝ^3

julia> op*ones(3,4)
3-element Vector{Float64}:
 4.0
 4.0
 4.0
	

MyLinOp(domain_type::Type, dim_in::Tuple, [domain_type::Type,] dim_out::Tuple, Fwd!::Function, Adj!::Function)

Construct a user defined LinearOperator by specifing its linear mapping Fwd! and its adjoint Adj!. The functions Fwd! and Adj must be in-place functions consistent with the given dimensions dim_in and dim_out and the domain and codomain types.

julia> n,m = 5,4;

julia> A = randn(n,m);

julia> op = MyLinOp(Float64, (m,),(n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))
A  ℝ^4 -> ℝ^5

julia> op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))
A  ℝ^4 -> ℝ^5
	

FiniteDiff([domain_type=Float64::Type,] dim_in::Tuple, direction = 1)
FiniteDiff(x::AbstractArray, direction = 1)

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns the discretized gradient over the specified direction obtained using forward finite differences.

julia> FiniteDiff(Float64,(3,))
δx  ℝ^3 -> ℝ^2

julia> FiniteDiff((3,4),2)
δy  ℝ^(3, 4) -> ℝ^(3, 3)

julia> all(FiniteDiff(ones(2,2,2,3),1)*ones(2,2,2,3) .== 0)
true
	

Variation([domaintype=Float64::Type,] dimin::Tuple) Variation(dims...) Variation(x::AbstractArray)

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns a matrix with its ith column consisting of the vectorized discretized gradient over the ith `direction obtained using forward finite differences.

julia> Variation(Float64,(10,2))
Ʋ  ℝ^(10, 2) -> ℝ^(20, 2)

julia> Variation(2,2,2)
Ʋ  ℝ^(2, 2, 2) -> ℝ^(8, 3)

julia> Variation(ones(2,2))*[1. 2.; 1. 2.]
4×2 Matrix{Float64}:
 0.0  1.0
 0.0  1.0
 0.0  1.0
 0.0  1.0

LBFGS(domain_type::Type,dim_in::Tuple, M::Integer)
LBFGS(dim_in::Tuple, M::Integer)
LBFGS(x::AbstractArray, M::Integer)

Construct a Limited-Memory BFGS LinearOperator with memory M. The memory of LBFGS can be updated using the function update!, where the current iteration variable and gradient (x, grad) and the previous ones (x_prev and grad_prev) are needed:

julia> L = LBFGS(Float64,(4,),5)
LBFGS  ℝ^4 -> ℝ^4

julia> update!(L,x,x_prev,grad,grad_prev); # update memory

julia> d = L*grad; # compute new direction

Use reset!(L) to cancel the memory of the operator.

reset!(L::LBFGS)

Cancels the memory of L.

update!(L::LBFGS, x, x_prex, grad, grad_prev)

See the documentation for LBFGS.

Exp([domain_type=Float64::Type,] dim_in::Tuple)

Creates the exponential non-linear operator with input dimensions dim_in:

\[e^{ \mathbf{x} }.\]

Pow([domain_type=Float64::Type,] dim_in::Tuple)

Elementwise power p non-linear operator with input dimensions dim_in.

Cos([domain_type=Float64::Type,] dim_in::Tuple)

Creates a cosine non-linear operator with input dimensions dim_in:

\[\cos (\mathbf{x} ).\]

Sin([domain_type=Float64::Type,] dim_in::Tuple)

Creates a sinusoid non-linear operator with input dimensions dim_in:

\[\sin( \mathbf{x} ).\]

Atan([domain_type=Float64::Type,] dim_in::Tuple)

Creates an inverse tangent non-linear operator with input dimensions dim_in:

\[\text{atan} ( \mathbf{x} ).\]

Tanh([domain_type=Float64::Type,] dim_in::Tuple)

Creates an hyperbolic tangent non-linear operator with input dimensions dim_in:

\[\text{tanh} ( \mathbf{x} ).\]

Sech([domain_type=Float64::Type,] dim_in::Tuple)

Creates an hyperbolic secant non-linear operator with input dimensions dim_in:

\[\text{sech} ( \mathbf{x} ).\]

Sigmoid([domain_type=Float64::Type,] dim_in::Tuple, γ = 1.)

Creates the sigmoid non-linear operator with input dimensions dim_in.

\[\sigma(\mathbf{x}) = \frac{1}{1+e^{-\gamma \mathbf{x} } }\]

SoftPlus([domain_type=Float64::Type,] dim_in::Tuple)

Creates the softplus non-linear operator with input dimensions dim_in.

\[\sigma(\mathbf{x}) = \log (1 + e^{x} )\]

SoftMax([domain_type=Float64::Type,] dim_in::Tuple)

Creates the softmax non-linear operator with input dimensions dim_in.

\[\sigma(\mathbf{x}) = \frac{e^{\mathbf{x} }}{ \sum e^{\mathbf{x} } }\]

DFT([domaintype=Float64::Type,] dimin::Tuple [,dims]; [normalization, flags, timelimit, numthreads]) DFT(dimin...; [normalization, flags, timelimit, numthreads]) DFT(x::AbstractArray [,dims]; [normalization, flags, timelimit, numthreads])

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns the N-dimensional Discrete Fourier Transform over dimensions dims of x.

Arguments:

• domain_type: The type of the input array. Defaults to Float64.
• dim_in: The dimensions of the input array. If dim_in is a tuple, it specifies the size of the input array. If dim_in is a type, it specifies the type of the input array.
• x: An input array. If provided, the dimensions of x will be used as dim_in.
• dims: The dimensions over which to perform the Fourier Transform. Defaults to all dimensions of x.
• flags: FFTW flags for the plan. Defaults to FFTW.ESTIMATE.
• normalization: The normalization scheme to use. The options are:
  • UNNORMALIZED: No normalization is applied (default).
  • ORTHO: Orthogonal normalization, scaling by 1/sqrt(N) both forward and backward transforms.
  • FORWARD: Forward normalization, scaling by 1/N.
  • BACKWARD: Backward normalization, scaling by 1/N.
• timelimit: The maximum time in seconds to spend on planning the FFT. Defaults to Inf, which means no time limit. If set to a finite value, the plan will be created within that time limit, potentially resulting in a less optimal plan.
• num_threads: The number of threads to use for the FFT. Defaults to the number of threads available in Julia. It should be set to 1 for single-threaded execution (e.g., when using within a @threads block).

julia> using FFTW, FFTWOperators

julia> DFT(Complex{Float64},(10,10))
ℱ  ℂ^(10, 10) -> ℂ^(10, 10)

julia> DFT(10,10)
ℱ  ℝ^(10, 10) -> ℂ^(10, 10)

julia> op = DFT(ones(3))
ℱ  ℝ^3 -> ℂ^3

julia> op*ones(3) ≈ FFTW.fft(ones(3))
true

IDFT([domaintype=Float64::Type,] dimin::Tuple [,dims]; [flags, timelimit, numthreads]) IDFT(dimin...; [flags, timelimit, numthreads]) IDFT(x::AbstractArray [,dims]; [flags, timelimit, numthreads])

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns the N-dimensional Inverse Discrete Fourier Transform over dimensions dims of x.

Arguments:

• domain_type: The type of the input array. Defaults to Float64.
• dim_in: The dimensions of the input array. If dim_in is a tuple, it specifies the size of the input array. If dim_in is a type, it specifies the type of the input array.
• x: An input array. If provided, the dimensions of x will be used as dim_in.
• dims: The dimensions over which to perform the Inverse Fourier Transform. Defaults to all dimensions of x.
• normalization: The normalization scheme to use. The options are:
  • UNNORMALIZED: No normalization is applied.
  • ORTHO: Orthogonal normalization, scaling by 1/sqrt(N) both forward and backward transforms.
  • FORWARD: Forward normalization, scaling by 1/N.
  • BACKWARD: Backward normalization, scaling by 1/N (default).
• flags: FFTW flags for the plan. Defaults to FFTW.ESTIMATE.
• timelimit: The maximum time in seconds to spend on planning the FFT. Defaults to Inf, which means no time limit. If set to a finite value, the plan will be created within that time limit, potentially resulting in a less optimal plan.
• num_threads: The number of threads to use for the FFT. Defaults to the number of threads available in Julia. It should be set to 1 for single-threaded execution (e.g., when using within a @threads block).

julia> using FFTW, FFTWOperators

julia> IDFT(Complex{Float64},(10,10))
ℱ⁻¹  ℂ^(10, 10) -> ℂ^(10, 10)

julia> IDFT(10,10)
ℱ⁻¹  ℂ^(10, 10) -> ℂ^(10, 10)

julia> op = IDFT(ones(ComplexF64, 3))
ℱ⁻¹  ℂ^3 -> ℂ^3

julia> op*ones(ComplexF64, 3) ≈ FFTW.ifft(ones(ComplexF64, 3))
true

RDFT([domain_type=Float64::Type,] dim_in::Tuple [,dims=1])
RDFT(dim_in...)
RDFT(x::AbstractArray [,dims=1])

Creates a LinearOperator which, when multiplied with a real array x, returns the DFT over the dimension dims, exploiting Hermitian symmetry.

julia> using FFTWOperators

julia> RDFT(Float64,(10,10))
ℱ  ℝ^(10, 10) -> ℂ^(6, 10)

julia> RDFT((10,10,10),2)
ℱ  ℝ^(10, 10, 10) -> ℂ^(10, 6, 10)
	

IRDFT([domain_type=Float64::Type,] dim_in::Tuple, d::Int, [,dims=1])
IRDFT(x::AbstractArray, d::Int, [,dims=1])

Creates a LinearOperator which, when multiplied with a complex array x, returns the IDFT over the dimension dims, exploiting Hermitian symmetry. Like in the function BASE.irfft, d must satisfy div(d,2)+1 == size(x,dims).

julia> using FFTWOperators

julia> A = IRDFT(Complex{Float64},(10,),19)
ℱ⁻¹  ℂ^10 -> ℝ^19

julia> A = IRDFT((5,10,8),19,2)
ℱ⁻¹  ℂ^(5, 10, 8) -> ℝ^(5, 19, 8)
	

DCT([domain_type=Float64::Type,] dim_in::Tuple)
DCT(dim_in...)
DCT(x::AbstractArray)

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns the N-dimensional Inverse Discrete Cosine Transform of x.

julia> using FFTWOperators

julia> DCT(Complex{Float64},(10,10))
ℱc  ℂ^(10, 10) -> ℂ^(10, 10)

julia> DCT(10,10)
ℱc  ℝ^(10, 10) -> ℝ^(10, 10)

julia> A = DCT(ones(3))
ℱc  ℝ^3 -> ℝ^3

julia> A*ones(3)
3-element Vector{Float64}:
 1.7320508075688772
 0.0
 0.0
	

IDCT([domain_type=Float64::Type,] dim_in::Tuple)
IDCT(dim_in...)
IDCT(x::AbstractArray)

Creates a LinearOperator which, when multiplied with an array x::AbstractArray{N}, returns the N-dimensional inverse Discrete Cosine Transform of x.

julia> using FFTWOperators

julia> IDCT(Complex{Float64},(10,10))
ℱc⁻¹  ℂ^(10, 10) -> ℂ^(10, 10)

julia> IDCT(10,10)
ℱc⁻¹  ℝ^(10, 10) -> ℝ^(10, 10)

julia> A = IDCT(ones(3))
ℱc⁻¹  ℝ^3 -> ℝ^3

julia> A*[1.;0.;0.]
3-element Vector{Float64}:
 0.5773502691896258
 0.5773502691896258
 0.5773502691896258

FFTShift([T::Type=Float64,] dim_in::Tuple, [dirs])
FFTShift(dim_in...)

Creates a LinearOperator that permutes the array like FFTW.fftshift over the given dirs. dirs must contain at least one dimension; each must be within 1:length(dim_in).

julia> using FFTWOperators

julia> A = FFTShift((4,)); x = collect(Float64, 1:4);

julia> A * x
4-element Vector{Float64}:
 3.0
 4.0
 1.0
 2.0

julia> B = IFFTShift((4,)); B * (A * x) == x
true

julia> A2 = FFTShift((2,2), (1,2)); A2 * reshape(1.0:4.0, 2, 2)
2×2 Matrix{Float64}:
 4.0  2.0
 3.0  1.0

IFFTShift([T::Type=Float64,] dim_in::Tuple, [dirs])
IFFTShift(dim_in...)

Creates a LinearOperator that permutes the array like FFTW.ifftshift over dirs. dirs must contain at least one dimension; each must be within 1:length(dim_in).

julia> using FFTWOperators

julia> A = FFTShift((4,)); B = IFFTShift((4,)); x = collect(1.0:4.0);

julia> B * (A * x) == x
true

julia> B2 = IFFTShift((2,2), (1,2)); B2 * [4.0 3.0; 2.0 1.0] == [1.0 2.0; 3.0 4.0]
true

SignAlternation([T::Type=Float64,] dim_in::Tuple, dirs)

Creates a LinearOperator that multiplies entries by -1 on indices where the parity sum across dirs is odd. dirs must contain at least one dimension; each must be within 1:length(dim_in).

Why is it useful?

Due to the properties of the discrete Fourier transform, when there is an even number of points along a dimension, alternating the sign of the entries along the domain of the Fourier transform (i.e., multiplying by -1 at every other index) is equivalent to a half-sample shift in the frequency domain. One can see this by applying the shift theorem of the Fourier transform: if $\mathcal{F}({x_n})_k = X_k$ then $\mathcal{F}({x_n \cdot e^{\frac{i 2\pi}{N}n m}})_k = \mathcal{F}({x_n \cdot (-1)^n })_k = X_{k - N/2}$ where N is the number of samples along that dimension.

Examples

julia> using FFTWOperators

julia> S = SignAlternation((4,), 1); S * collect(1.0:4.0)
4-element Vector{Float64}:
  1.0
 -2.0
  3.0
 -4.0

julia> S2 = SignAlternation((2,2), (1,2)); S2 * ones(2,2)
2×2 Matrix{Float64}:
  1.0  -1.0
 -1.0   1.0

fftshift_op(op::AbstractOperator; domain_shifts::Tuple=(), codomain_shifts::Tuple=())

Applies FFTShift to the domain and/or codomain of op, depending on domain_shifts and codomain_shifts. If the shifted dimensions all have even length and op is a (possibly modified) DFT/IDFT, the FFTShift is replaced by a SignAlternation on the other side of the DFT/IDFT.

Examples

julia> using FFTWOperators, FFTW

julia> x = rand(15);

julia> F = fftshift_op(DFT(15); domain_shifts=(1,))
ℱ*⇉  ℝ^15 -> ℂ^15

julia> F * x ≈ FFTW.fft(FFTW.fftshift(x))
true

julia> F = fftshift_op(DFT(15); codomain_shifts=(1,))
⇉*ℱ  ℝ^15 -> ℂ^15

julia> F * x ≈ FFTW.fftshift(FFTW.fft(x))
true

julia> F = fftshift_op(DFT(15); domain_shifts=(1,), codomain_shifts=(1,))
Π  ℝ^15 -> ℂ^15

julia> F * x ≈ FFTW.fftshift(FFTW.fft(FFTW.fftshift(x)))
true

julia> F = fftshift_op(DFT(16); codomain_shifts=(1,)) # note that 16 is even, so we get a SignAlternation (±)
ℱ*±  ℝ^16 -> ℂ^16

ifftshift_op(op::AbstractOperator; domain_shifts::Tuple=(), codomain_shifts::Tuple=())

Applies IFFTShift to the domain and/or codomain of op, depending on domain_shifts and codomain_shifts. If the shifted dimensions all have even length and op is a (possibly modified) DFT/IDFT, the IFFTShift is replaced by a SignAlternation on the other side of the DFT/IDFT.

Examples

julia> using FFTWOperators, FFTW

julia> x = rand(ComplexF64, 15);

julia> F = ifftshift_op(IDFT(15); domain_shifts=(1,))
ℱ⁻¹*⇇  ℂ^15 -> ℂ^15

julia> F * x ≈ FFTW.ifft(FFTW.ifftshift(x))
true

julia> F = ifftshift_op(IDFT(15); codomain_shifts=(1,))
⇇*ℱ⁻¹  ℂ^15 -> ℂ^15

julia> F * x ≈ FFTW.ifftshift(FFTW.ifft(x))
true

julia> F = ifftshift_op(IDFT(15); domain_shifts=(1,), codomain_shifts=(1,))
Π  ℂ^15 -> ℂ^15

julia> F * x ≈ FFTW.ifftshift(FFTW.ifft(FFTW.ifftshift(x)))
true

julia> F = ifftshift_op(IDFT(16); codomain_shifts=(1,)) # note that 16 is even, so we get a SignAlternation (±)
ℱ⁻¹*±  ℂ^16 -> ℂ^16

alternate_sign(x[, dirs...])

Returns a copy with sign alternation across specified dimensions.

julia> using FFTWOperators

julia> alternate_sign(collect(1.0:4.0), 1)
4-element Vector{Float64}:
  1.0
 -2.0
  3.0
 -4.0

    alternate_sign!(x[, dirs...])

In-place sign alternation across specified dimensions. Flips the sign where the parity sum across dirs is odd. The provided dirs must be within 1:ndims(x), and sorted in ascending order. They can also be provided as a tuple.

julia> using FFTWOperators

julia> v = collect(1.0:4.0);

julia> alternate_sign!(v, 1)
4-element Vector{Float64}:
  1.0
 -2.0
  3.0
 -4.0

julia> M = ones(2,2); alternate_sign!(M, 1, 2)
2×2 Matrix{Float64}:
  1.0  -1.0
 -1.0   1.0

julia> M = ones(2,2); alternate_sign!(M, (1, 2))
2×2 Matrix{Float64}:
  1.0  -1.0
 -1.0   1.0

alternate_sign!(y, x[, dirs...])

Out-of-place variant: apply sign alternation of x across dirs and store the result in y. y and x must have the same size. The provided dirs must be within 1:ndims(x), and sorted in ascending order. They can also be provided as a tuple.

julia> using FFTWOperators

julia> x = reshape(1.0:4.0, 2, 2); y = similar(x);

julia> alternate_sign!(y, x, 1, 2)
2×2 Matrix{Float64}:
  1.0  -3.0
 -2.0   4.0

julia> alternate_sign!(y, x, (1, 2))
2×2 Matrix{Float64}:
  1.0  -3.0
 -2.0   4.0

Conv([domain_type=Float64::Type,] dim_in::Tuple, h::AbstractVector)
Conv(x::AbstractVector, h::AbstractVector)

Creates a LinearOperator which, when multiplied with an array x::AbstractVector, returns the convolution between x and h. Uses conv and hence FFT algorithm.

julia> using DSPOperators

julia> Conv((10,),randn(5))
★  ℝ^10 -> ℝ^14
	

Xcorr([domain_type=Float64::Type,] dim_in::Tuple, h::AbstractVector)
Xcorr(x::AbstractVector, h::AbstractVector)

Creates a LinearOperator which, when multiplied with an array x::AbstractVector, returns the cross correlation between x and h. Uses xcorr from DSP.jl.

Examples

julia> using DSPOperators

julia> Xcorr(Float64, (10,), [1.0, 0.5, 0.2])
◎  ℝ^10 -> ℝ^19

Filt([domain_type=Float64::Type,] dim_in::Tuple, b::AbstractVector, [a::AbstractVector,])
Filt(x::AbstractVector, b::AbstractVector, [a::AbstractVector,])

Creates a LinearOperator which, when multiplied with an array x::AbstractVector, returns a vector y filtered by an IIR filter of coefficients b and a. If only b is provided a FIR is used to comute y instead.

julia> using DSPOperators

julia> Filt(Float64, (10,), [1.0, 2.0, 1.0], [1.0, -1.0])
IIR  ℝ^10 -> ℝ^10
	

MIMOFilt([domain_type=Float64::Type,] dim_in::Tuple, B::Vector{AbstractVector}, [A::Vector{AbstractVector},])

MIMOFilt(x::AbstractMatrix, b::Vector{AbstractVector}, [a::Vector{AbstractVector},])

Creates a LinearOperator which, when multiplied with a matrix X, returns a matrix Y. Here a Multiple Input Multiple Output system is evaluated: the columns of X and Y represent the input signals and output signals respectively.

\[\mathbf{y}_i = \sum_{j = 1}^{M} \mathbf{h}_{i,j} * \mathbf{x}_j\]

where $\mathbf{y}_i$ and $\mathbf{x}_j$ are the $i$-th and $j$-th columns of the output Y and input X matrices respectively.

The filters $\mathbf{h}_{i,j}$ can be represented either by providing coefficients B and A (IIR) or B alone (FIR). These coefficients must be given in a Vector of Vectors.

For example for a 3 by 2 MIMO system (i.e. size(X,2) == 3 inputs and size(Y,2) == 2 outputs) B must be:

B = [b11, b12, b13, b21, b22, b23]

where bij are vector containing the filter coeffients of h_{i,j}.

julia> using DSPOperators

julia> m,n = 10,3; #time samples, number of inputs

julia> B  = [[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.], ];

julia> #B = [   b11   ,     b12   ,    b13   ,   b21    ,   b22,       b23    , ]


julia> A  = [[1.;1.;1.],[2.;2.;2.],[      3.],[      4.],[      5.],[      6.], ];

julia> #A = [   a11   ,     a12   ,    a13   ,   a21    ,   a22,       a23    , ]


julia> op = MIMOFilt(Float64, (m,n), B, A)
※  ℝ^(10, 3) -> ℝ^(10, 2)

julia> X = randn(m,n); #input signals

julia> Y = op*X;       #output signals

julia> using DSP

julia> Y[:,1] ≈ filt(B[1],A[1],X[:,1])+filt(B[2],A[2],X[:,2])+filt(B[3],A[3],X[:,3])
true

julia> Y[:,2] ≈ filt(B[4],A[4],X[:,1])+filt(B[5],A[5],X[:,2])+filt(B[6],A[6],X[:,3])
true
	

ZeroPad([domain_type::Type,] dim_in::Tuple, zp::Tuple)
ZeroPad(x::AbstractArray, zp::Tuple)

Create a LinearOperator which, when multiplied to an array x of size dim_in, returns an expanded array y of size dim_in .+ zp where y[1:dim_in[1], 1:dim_in[2] ... ] = x and zero elsewhere.

julia> Z = ZeroPad((2,2),(0,2))
[I;0]  ℝ^(2, 2) -> ℝ^(2, 4)

julia> Z*ones(2,2)
2×4 Matrix{Float64}:
 1.0  1.0  0.0  0.0
 1.0  1.0  0.0  0.0
	

NFFTOp(image_size::NTuple{D,Int}, trajectory::AbstractArray{T}, dcf::AbstractArray; threaded::Bool=true, kwargs...)

Create a non-uniform fast Fourier transform operator [1]. The operator is created with a given image size, trajectory, and density compensation function (dcf). The dcf is used to correct for the non-uniform sample density of the trajectory. The operator can be used to transform images to k-space and back.

<em>To use the operator, the NFFT package must be explicitly imported!</em>

Arguments

• image_size::NTuple{D,Int}: The size of the image to transform.
• trajectory::AbstractArray{T}: The trajectory of the samples in k-space. The first dimension of the trajectory must match the number of image dimensions. The trajectory must have at least two dimensions.
• dcf::AbstractArray: The density compensation function. The shape of the trajectory from the second dimension must match the shape of the dcf array. The element type of the trajectory must match the element type of the dcf array. This argument is optional and defaults to nothing. If nothing is passed, the dcf will be estimated using the sample density compensation method [2].
• threaded::Bool=true: Whether to use threading when applying the operator. Defaults to true.
• dcf_estimation_iterations::Union{Nothing,Int}=nothing: The number of iterations to use when estimating the dcf. Defaults to 20. This argument is only used if dcf is not provided.
• dcf_correction_function::Function=identity: A correction function to apply to the estimated dcf. Defaults to the identity function. This argument is only used if dcf is not provided.
• kwargs...: Additional keyword arguments to pass to the NFFTPlan constructor.

References

1. Fessler, J. A., & Sutton, B. P. (2003). Nonuniform fast Fourier transforms using min-max interpolation.

IEEE Transactions on Signal Processing, 51(2), 560-574.

1. Pipe, J. G., & Menon, P. (1999). Sampling density compensation in MRI: rationale and an iterative numerical solution.

Examples

julia> using NFFTOperators

julia> image_size = (128, 128);

julia> trajectory = rand(2, 128, 50) .- 0.5;

julia> dcf = rand(128, 50);

julia> op = NFFTOp(image_size, trajectory, dcf)
𝒩  ℂ^(128, 128) -> ℂ^(128, 50)

julia> image = rand(ComplexF64, image_size);

julia> ksp = op * image;

julia> image_reconstructed = op' * ksp;

WaveletOp(wavelet::DiscreteWavelet, dim_in::Integer)
WaveletOp(wavelet::DiscreteWavelet, dim_in::Tuple)

Creates a LinearOperator which, when multiplied with a vector x::AbstractVector, returns the wavelet transform of x using the given wavelet and levels.

julia> using WaveletOperators

julia> W = WaveletOp(wavelet(WT.db4), 4)
𝒲  ℝ^4 -> ℝ^4

julia> W * ones(4)
4-element Vector{Float64}:
  2.0
 -5.551115123125783e-17
 -8.326672684688674e-17
 -8.326672684688674e-17