function AbstractNFFTs.deconvolve!(p::NFFTPlan{T,D,R}, f::AbstractArray{U,D},
                      g::StridedArray{Complex{T},D}) where {T,D,R,U}
  dstart = Val(p.dims[1])
  dend = Val(p.dims[end])
  return _deconvolve!(p, f, g, dstart, dend)
end

@generated function _deconvolve!(p::NFFTPlan{T,D,R}, f::AbstractArray{U,D},
                                  g::StridedArray{Complex{T},D}, dstart::Val{DS}, 
                                  dend::Val{DE}) where {T,D,R,U,DS,DE}
    quote
        @nloops $D l f d->begin
            # preexpr
            if $DS ≤ d ≤ $DE
                offset_d = round(Int, p.n[d] - p.N[d]÷2) - 1
                gidx_d = rem(l_d+offset_d, p.n[d]) + 1
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            v = (@nref $D f l)
            @nexprs $(DE-DS+1) d -> v *= p.windowHatInvLUT[d][l_{d+$DS-1}]
            (@nref $D g gidx) = v
        end
    end
end

function AbstractNFFTs.deconvolve_transpose!(p::NFFTPlan{T,D,R},
                 g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {T,D,R,U}
  dstart = Val(p.dims[1])
  dend = Val(p.dims[end])
  return _deconvolve_transpose!(p, g, f, dstart, dend)
end

@generated function _deconvolve_transpose!(p::NFFTPlan{T,D,R},
           g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}, dstart::Val{DS}, 
           dend::Val{DE}) where {T,D,R,U,DS,DE}
    quote
        @nloops $D l f d->begin
            # preexpr
            if $DS ≤ d ≤ $DE
                offset_d = round(Int, p.n[d] - p.N[d]÷2) - 1
                gidx_d = rem(l_d+offset_d, p.n[d]) + 1
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            v =  (@nref $D g gidx)
            @nexprs $(DE-DS+1) d -> v *= p.windowHatInvLUT[d][l_{d+$DS-1}]
            (@nref $D f l) = v
        end
    end
end

function AbstractNFFTs.convolve!(p::NFFTPlan{T,D,R}, g::AbstractArray{Complex{T},D},
      fHat::StridedArray{U,R}) where {T,D,R,U}
  l1 = Val(p.dims[1]-1)
  l2 = Val(length(p.dims))
  l3 = Val(D-p.dims[end])
  return _convolve!(p, g, fHat, l1, l2, l3)
end

@generated function _convolve!(p::NFFTPlan{T,D,R}, g::AbstractArray{Complex{T},D},
                   fHat::StridedArray{U,R}, l1::Val{L1}, l2::Val{L2}, l3::Val{L3}) where {T,D,R,U,L1,L2,L3}
  quote
    fill!(fHat, zero(T))
    scale = Int(p.params.LUTSize/(p.params.m+2))

    for k in 1:p.M
      ### The first L1 loops contain no NFFT
      @nloops_ $L1 d->l_{d} d->1:size(g,d) d->begin
         # preexpr
         gidx_d = l_d
         fidx_d = l_d
      end begin
          # bodyexpr
          @nexprs 1 d -> fidx_{$L1+1} = k
          @nexprs 1 d -> prodWin_{$L2} = one(T)
          @nexprs $L2 d -> xscale_d = p.x[d,k] * p.n[d+$L1]
          @nexprs $L2 d -> c_d = floor(Int, xscale_d)
          ### The next L2 loops contain the NFFT
          @nloops_ $L2 d->l_{d+$L1} d->(c_d-p.params.m):(c_d+p.params.m) d->begin
            # preexpr
            gidx_{d+$L1} = rem(l_{d+$L1}+p.n[d+$L1], p.n[d+$L1]) + 1
            idx = abs((xscale_d - l_{d+$L1} )*scale) + 1
            idxL = floor(idx)
            idxInt = Int(idxL)
            prodWin_{d-1} = prodWin_d * (p.windowLinInterp[idxInt] + ( idx-idxL ) * (p.windowLinInterp[idxInt+1] - p.windowLinInterp[idxInt]))
          end begin
             # bodyexpr
             ### The last L3 loops contain no NFFT
             @nloops_ $L3 d->l_{d+$L1+$L2} d->1:size(g,d+$L1+$L2) d->begin
                 # preexpr
                 gidx_{d+$L1+$L2} = l_{d+$L1+$L2}
                 fidx_{d+$L1+1} = l_{d+$L1+$L2}
              end begin  
                (@nref $R fHat fidx) += (@nref $D g gidx) * prodWin_0
              end 
          end
      end
    end
  end
end

function AbstractNFFTs.convolve_transpose!(p::NFFTPlan{T,D,R}, fHat::AbstractArray{U,R},
        g::StridedArray{Complex{T},D}) where {T,D,R,U}
  l1 = Val(p.dims[1]-1)
  l2 = Val(length(p.dims))
  l3 = Val(D-p.dims[end])
  return _convolve_transpose!(p, fHat, g, l1, l2, l3)
end

@generated function _convolve_transpose!(p::NFFTPlan{T,D,R}, fHat::AbstractArray{U,R},
            g::StridedArray{Complex{T},D}, l1::Val{L1}, l2::Val{L2}, l3::Val{L3}) where {T,D,R,U,L1,L2,L3}
  quote
    fill!(g, zero(T))
    scale = Int(p.params.LUTSize/(p.params.m+2))

    for k in 1:p.M
      ### The first L1 loops contain no NFFT
      @nloops_ $L1 d->l_{d} d->1:size(g,d) d->begin
         # preexpr
         gidx_d = l_d
         fidx_d = l_d
      end begin
          # bodyexpr
          @nexprs 1 d -> fidx_{$L1+1} = k
          @nexprs 1 d -> prodWin_{$L2} = one(T)
          @nexprs $L2 d -> xscale_d = p.x[d,k] * p.n[d+$L1]
          @nexprs $L2 d -> c_d = floor(Int, xscale_d)
          ### The next L2 loops contain the NFFT
          @nloops_ $L2 d->l_{d+$L1} d->(c_d-p.params.m):(c_d+p.params.m) d->begin
            # preexpr
            gidx_{d+$L1} = rem(l_{d+$L1}+p.n[d+$L1], p.n[d+$L1]) + 1
            idx = abs((xscale_d - l_{d+$L1})*scale) + 1
            idxL = floor(idx)
            idxInt = Int(idxL)
            prodWin_{d-1} = prodWin_d * (p.windowLinInterp[idxInt] + ( idx-idxL ) * (p.windowLinInterp[idxInt+1] - p.windowLinInterp[idxInt]))
          end begin
             # bodyexpr
             ### The last L3 loops contain no NFFT
             @nloops_ $L3 d->l_{d+$L1+$L2} d->1:size(g,d+$L1+$L2) d->begin
                 # preexpr
                 gidx_{d+$L1+$L2} = l_{d+$L1+$L2}
                 fidx_{d+$L1+1} = l_{d+$L1+$L2}
              end begin  
                (@nref $D g gidx) += (@nref $R fHat fidx) * prodWin_0
              end 
          end
      end
    end
  end
end