/
RDFT.jl
92 lines (69 loc) · 2.08 KB
/
RDFT.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
export RDFT
"""
`RDFT([domainType=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
julia> RDFT(Float64,(10,10))
ℱ ℝ^(10, 10) -> ℂ^(6, 10)
julia> RDFT((10,10,10),2)
ℱ ℝ^(10, 10, 10) -> ℂ^(10, 6, 10)
```
"""
struct RDFT{T <:Number,
N,
T1<:AbstractFFTs.Plan,
T2<:AbstractFFTs.Plan,
T3<:AbstractArray{Complex{T},N}
} <: LinearOperator
dim_in::NTuple{N,Int}
dim_out::NTuple{N,Int}
A::T1
At::T2
b2::T3
y2::T3
Zp::ZeroPad{Complex{T},N}
end
# Constructors
#standard constructor
function RDFT(x::AbstractArray{T,N}, dims::Int=1) where {T<:Real,N}
A = plan_rfft(x,dims)
b2 = zeros(complex(T),size(x))
y2 = zeros(complex(T),size(x))
At = plan_bfft(y2,dims)
dim_in = size(x)
dim_out = ()
for i = 1:N
dim_out = i == dims ? (dim_out...,div(dim_in[i],2)+1) : (dim_out...,dim_in[i])
end
Z = ZeroPad(Complex{T},dim_out,size(b2))
RDFT{T,N,typeof(A),typeof(At),typeof(b2)}(dim_in,dim_out,A,At,b2,y2,Z)
end
RDFT(T::Type,dim_in::NTuple{N,Int},dims::Int=1) where {N} = RDFT(zeros(T,dim_in),dims)
RDFT(dim_in::NTuple{N,Int},dims::Int=1) where {N} = RDFT(zeros(dim_in),dims)
RDFT(dim_in::Vararg{Int}) = RDFT(dim_in)
RDFT(T::Type,dim_in::Vararg{Int}) = RDFT(T,dim_in)
# Mappings
function mul!(y::T3,
L::RDFT{T,N,T1,T2,T3},
b::T4) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}}
mul!(y,L.A,b)
end
function mul!(y::T4,
L::AdjointOperator{RDFT{T,N,T1,T2,T3}},
b::T3) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}}
A = L.A
mul!(A.b2,A.Zp,b)
mul!(A.y2,A.At,A.b2)
y .= real.(A.y2)
return y
end
# Properties
size(L::RDFT) = (L.dim_out,L.dim_in)
fun_name(A::RDFT) = "ℱ"
domainType(L::RDFT{T,N}) where {T,N} = T
codomainType(L::RDFT{T,N}) where {T,N} = Complex{T}
is_AAc_diagonal( L::RDFT) = false #TODO but might be true?
is_invertible( L::RDFT) = true
is_full_row_rank( L::RDFT) = true