"Kernel density estimation" is a convolution of what? I am trying to get a better understanding of kernel density estimation.
Using the definition from Wikipedia: https://en.wikipedia.org/wiki/Kernel_density_estimation#Definition
$
 \hat{f_h}(x) = \frac{1}{n}\sum_{i=1}^n K_h (x - x_i) \quad = \frac{1}{nh} \sum_{i=1}^n K\Big(\frac{x-x_i}{h}\Big)
$
Let's take $K()$ to be a rectangular function which gives $1$ if $x$ is between $-0.5$ and $0.5$ and $0$ otherwise, and $h$ (window size) to be 1.
I understand that the density is a convolution of two functions, but I am not sure I know how to define these two functions. One of them should (probably) be a function of the data which, for every point in R, tells us how many data points we have in that location (mostly $0$). And the other function should probably be some modification of the kernel function, combined with the window size. But I am not sure how to define it. 
Any suggestions?
Bellow is an example R code which (I suspect) replicates the settings I defined above (with a mixture of two Gaussians and $n=100$), on which I hope to see a "proof" that the functions to be convoluted are as we suspect.
# example code:
set.seed(2346639)
x <- c(rnorm(50), rnorm(50,2))
plot(density(x, kernel='rectangular', width=1, n = 10**4))
rug(x)


 A: Corresponding to any batch of data $X = (x_1, x_2, \ldots, x_n)$ is its "empirical density function"
$$f_X(x) = \frac{1}{n}\sum_{i=1}^{n} \delta(x-x_i).$$
Here, $\delta$ is a "generalized function."  Despite that name, it isn't a function at all: it's a new mathematical object that can be used only within integrals.  Its defining property is that for any function $g$ of compact support that is continuous in a neighborhood of $0$,
$$\int_{\mathbb{R}}\delta(x) g(x) dx = g(0).$$
(Names for $\delta$ include "atomic" or "point" measure and "Dirac delta function."  In the following calculation this concept is extended to include functions $g$ which are continuous from one side only.)
Justifying this characterization of $f_X$ is the observation that
$$\eqalign{
\int_{-\infty}^{x} f_X(y) dy 
&= \int_{-\infty}^{x} \frac{1}{n}\sum_{i=1}^{n} \delta(y-x_i)dy \\
&= \frac{1}{n}\sum_{i=1}^{n} \int_{-\infty}^{x}  \delta(y-x_i)dy \\
&= \frac{1}{n}\sum_{i=1}^{n} \int_{\mathbb{R}} I(y\le x) \delta(y-x_i)dy \\
&= \frac{1}{n}\sum_{i=1}^{n} I(x_i \le x) \\
&= F_X(x)
}$$
where $F_X$ is the usual empirical CDF and $I$ is the usual characteristic function (equal to $1$ where its argument is true and $0$ otherwise).  (I skip an elementary limiting argument needed to move from functions of compact support to functions defined over $\mathbb{R}$; because $I$ only needs to be defined for values within the range of $X$, which is compact, this is no problem.)
The convolution of $f_X(x)$ with any other function $k$ is given, by definition, as
$$\eqalign{
(f_X * k)(x) &= \int_{\mathbb{R}} f_X(x - y) k(y) dy \\
             &=\int_{\mathbb{R}}  \frac{1}{n}\sum_{i=1}^{n} \delta(x-y-x_i) k(y) dy \\
             &= \frac{1}{n}\sum_{i=1}^{n}\int_{\mathbb{R}}  \delta(x-y-x_i) k(y) dy \\
             &=\frac{1}{n}\sum_{i=1}^{n} k(x_i-x).
}$$
Letting $k(x) = K_h(-x)$ (which is the same as $K_h(x)$ for symmetric kernels--and most kernels are symmetric) we obtain the claimed result: the Wikipedia formula is a convolution.
