Let $a,b \in \mathbb{R}$ with $a < b$. Let also $P = \{x_1, \ldots, x_n\}$ (with $n > 1$) be a finite subset of $[a,b]$ with all distinct elements ($x_1 < x_2 < \cdots < x_n$). Let $C[a, b]$ be the real vector space of all continuous maps $f: [a, b] \to \mathbb{R}$. Of course $C[a, b]$ is also an inner product space with inner product $\langle f, g\rangle = \int_a^b f(x)g(x) dx$.
Is there an orthonormal basis for the space $C[a, b]$ whose (infinitely many) elements $\{e_k\}_{k \in \mathbb{N}}$ satisfy the constraint that
(*) for every $k \in \mathbb{N}$ there exists one index $i \in \{1,\ldots,n\}$ such that $e_k(x_j) = \delta_{ij}$ for all $j \in \{1,\ldots,n\}$ ?
Some polynomials exhibiting this property are the Lagrange polynomials $P_j(x) = \prod_{\substack{i=1 \\ i\neq j}}^{n} \frac{x-x_j}{x_i-x_j}$ (for $j = 1,\ldots, n$). However they are only finitely many so they cannot be a basis for the infinite-dimensional space $C[a,b]$.
The functions $P_1, \ldots, P_n$ belong to $C[a, b]$ and are linearly independent but not orthogonal. It should be possible to extend the set $\{P_1, \ldots, P_n\}$ to a base for $C[a, b]$ and then othonormalize it with Gram-Schmit. However extension and orthogonalization probably do not preserve the property (*).
MOTIVATION: When performing functional Principal Component Analysis in an infinite-dimensional function space one starts with a basis of the function space and finds a finite portion of that basis generating a finite-dimensional subspace with some desired properties. I would use a special basis with the property (*).
Best Answer
This is less difficult than you are making it. All you need to do is demarcate non-overlapping intervals $I_j$ about each of the $x_j$. Then for each $x_j$ define a function $e_j$ with the properties
There are many ways to accomplish this. For example,
For a given parameter $h$, define the function $g_h:\Bbb R \to \Bbb R$ by $$g_h(x) = \begin{cases}0& x\notin (-1,1)\\ 2h(x+1) & x\in(-1,-1/2]\\ 2(1-h)(x+1/2)+h & x \in(-1/2,0]\\ 2(h-1)x + 1 & x \in(0, 1/2]\\ -2h(x-1/2) + h & x \in (1/2, 1)\end{cases}$$Note that $g_h(0) = 1$ and $g_h(\pm1/2) = h$.
for each $1 \le j \le n$, select an $0 <\epsilon_j < \min \left\{\frac{x_j-x_{j-1}}2,\frac{x_{j+1}-x_j}2, 1\right\}$ (for $j=1,n$, dropping the term that is undefined). Set $I_j = (x_j -\epsilon_j, x_j+\epsilon_j) \cap [a,b]$.
For given parameter $h_k$, define $$e_k(x) = g_{h_k}\left(\dfrac{x - x_k}{\epsilon_k}\right), x \in [a,b]$$Note that $e_k(x_k) = 1$.
Choose each $h_k$ so that $\int_a^b e_k^2(x)\,dx = 1$. (Can you see why this will always be possible?)
Because the support (where they are non-zero) of the $e_k$ are each limited to the intervals $I_k$, which do not overlap, we have $\langle e_i, e_j\rangle = \delta_{ij}$ as desired.
The set $\{e_k \mid 1\le k \le n\}$ spans a subspace $V$ of $C[a,b]$. Choose any orthonormal basis $\{e_k\}_{k > n}$ for $V^\perp$. Then the combined basis $\{e_k\}_{k =1}^\infty$ is an orthonormal basis for $C[a,b]$ having all the properties you wanted.