diff --git a/src/bivariate.jl b/src/bivariate.jl
index d54c41ae..a6e18372 100644
--- a/src/bivariate.jl
+++ b/src/bivariate.jl
@@ -19,7 +19,7 @@ mutable struct BivariateKDE{Rx<:AbstractRange,Ry<:AbstractRange} <: AbstractKDE
     "Second coordinate of gridpoints for evaluating the density."
     y::Ry
     "Kernel density at corresponding gridpoints `Tuple.(x, permutedims(y))`."
-    density::Matrix{Float64}
+    density::AbstractMatrix{}
 end
 
 function kernel_dist(::Type{D},w::Tuple{Real,Real}) where D<:UnivariateDistribution
@@ -54,7 +54,7 @@ function tabulate(data::Tuple{RealVector, RealVector}, midpoints::Tuple{Rx, Ry},
     sx, sy = step(xmid), step(ymid)
 
     # Set up a grid for discretized data
-    grid = zeros(Float64, nx, ny)
+    grid = zeros(eltype(xdata),nx,ny)
     ainc = 1.0 / (sum(weights)*(sx*sy)^2)
 
     # weighted discretization (cf. Jones and Lotwick)
@@ -91,11 +91,12 @@ function conv(k::BivariateKDE, dist::Tuple{UnivariateDistribution,UnivariateDist
             ft[i+1,j+1] *= cf(distx,i*cx)*cf(disty,min(j,Ky-j)*cy)
         end
     end
-    dens = irfft(ft, Kx)
+    # i = 0:size(ft,1)-1
+    # j = 0:size(ft,2)-1
+    # ft = ft .* ( cf.(distx,i*cx) * cf.(disty,min.(j,Ky-j)*cy)' )
 
-    for i = 1:length(dens)
-        dens[i] = max(0.0,dens[i])
-    end
+    dens = irfft(ft, Kx)
+    dens = max.(0.0,dens)
 
     # Invert the Fourier transform to get the KDE
     BivariateKDE(k.x, k.y, dens)
diff --git a/src/univariate.jl b/src/univariate.jl
index ac92b1d8..a0132e73 100644
--- a/src/univariate.jl
+++ b/src/univariate.jl
@@ -17,7 +17,7 @@ mutable struct UnivariateKDE{R<:AbstractRange} <: AbstractKDE
     "Gridpoints for evaluating the density."
     x::R
     "Kernel density at corresponding gridpoints `x`."
-    density::Vector{Float64}
+    density::AbstractVector{}
 end
 
 # construct kernel from bandwidth
@@ -101,7 +101,7 @@ function tabulate(data::RealVector, midpoints::R, weights::Weights=default_weigh
     s = step(midpoints)
 
     # Set up a grid for discretized data
-    grid = zeros(Float64, npoints)
+    grid = zeros(eltype(data),npoints)
     ainc = 1.0 / (sum(weights)*s*s)
 
     # weighted discretization (cf. Jones and Lotwick)
@@ -125,22 +125,12 @@ function conv(k::UnivariateKDE, dist::UnivariateDistribution)
     K = length(k.density)
     ft = rfft(k.density)
 
-    # Convolve fft with characteristic function of kernel
-    # empirical cf
-    #  = \sum_{n=1}^N e^{i*t*X_n} / N
-    #  = \sum_{k=0}^K e^{i*t*(a+k*s)} N_k / N
-    #  = e^{i*t*a} \sum_{k=0}^K e^{-2pi*i*k*(-t*s*K/2pi)/K} N_k / N
-    #  = A * fft(N_k/N)[-t*s*K/2pi + 1]
     c = -twoπ/(step(k.x)*K)
-    for j = 0:length(ft)-1
-        ft[j+1] *= cf(dist,j*c)
-    end
+    j = 0:length(ft)-1
+    ft = ft .* cf.(dist,j*c)
 
     dens = irfft(ft, K)
-    # fix rounding error.
-    for i = 1:K
-        dens[i] = max(0.0,dens[i])
-    end
+    dens = max.(0.0,dens)
 
     # Invert the Fourier transform to get the KDE
     UnivariateKDE(k.x, dens)