Let $Q$ be a finite set with $|Q| = q$. Let $n \in \mathbb N$ be fixed, and let $a, b \in Q^n$ two elements with $d(a, b) = k$, where $d$ is the Hamming distance, i.e. the number of $1 ≤ i ≤ n$ with $a_i ≠ b_i$. We then define the set:
$$c_k(r, s) := \left| \left\{ c \in Q^n: d(a, c) = r \text{ and } d(b, c) = s\right\}\right|$$
I now want to show that:
$$P_r(i) P_s(i) = \sum_{k=0}^n c_k(r, s) P_k(i)$$
where $P_r(i), P_s(i)$ are the Krawtchuck polynomials, i.e. $P_r(i) = \sum_{j=0}^r (-1)^j (q - 1)^{r-j} \pmatrix{i \\ j} \pmatrix{n - i \\ r - j}$.
I must admit that I couldn't really get started so far. I've noticed that the set $c_k(r, s)$ is independent of the choice of $a, b$ as long as they satisfy the condition $d(a, b) = k$ (otherwise, it wouldn't make sense to not include $a, b$ as parameters for $c$).
I don't really know though how I can get from the Krawtchuck polynomials to that formula. I know some basic properties of the Krawtchuck polynomials like their orthogonality relation(s):
$\displaystyle \sum_{i=0}^n P_r(i) P_i(s) = q^n \delta_{r, s} \\ \sum_{i=0}^n \pmatrix{n \\ i} (q - 1)^i P_r(i) P_s(i) = q^n \pmatrix{n \\ r} (q - 1)^r \delta_{r, s} $
(where $\delta_{r, s}$ is the Kronecker-Delta, and where $q ≥ 2$). And I suspect that I somehow have to use them and utilize some smart properties or observations about these $c_k(r, s)$.
Without loss of generality, we may assume that $Q$ is the ring $\mathbb{Z}_q=\mathbb{Z}/q\mathbb{Z}$. We need the following lemma.
Lemma Let $\omega$ be a primitive $q$-th root of unity in $\mathbb{C}$ and let $x\in \mathbb{Z}_q^n$ with $wt(x)=i$, where $wt(x)$ denotes the Hamming weight of $x$. Then, we have $$\sum\limits_{y\in \mathbb{Z}_q^n\\wt(y)=k}\omega^{<x,y>}=P_k(i),\qquad\qquad (1)$$ where $\langle x,y\rangle$ denotes the usual inner product $\sum\limits_{l=1}^n x_l y_l$ of $x$ and $y$.
By using the lemma, we get for a fixed $x\in \mathbb{Z}_q^n$ with $ wt(x)=i$ that $$P_r(i) P_s(i) = \sum\limits_{y\in \mathbb{Z}_q^n\\wt(y)=r}\omega^{<x,y>} \sum\limits_{z\in \mathbb{Z}_q^n\\wt(z)=s}\omega^{<x,z>} = \sum\limits_{y\in \mathbb{Z}_q^n\\wt(y)=r} \sum\limits_{z\in \mathbb{Z}_q^n\\wt(z)=s}\omega^{<x,y+z>}. $$ Now, we use the numbers $c_k(r,s)$ (which are called intersection numbers). For all $u\in \mathbb{Z}_q^n$ with $wt(u)=d(u,0)=k$, we have $$c_k(r,s)=|\{ y\in\mathbb{Z}_q^n : wt(y)=d(y,0)=r \text{ and } d(y,u)=s \}|. $$ Therefore, $$\sum\limits_{z\in \mathbb{Z}_q^n\\wt(z)=s}\omega^{<x,y+z>}=\sum\limits_{k=0}^n c_k(r,s) \sum\limits_{u\in \mathbb{Z}_q^n\\wt(u)=k} \omega^{<x,u>}$$ and again by using the lemma, we get the desired equation $$P_r(i) P_s(i) = \sum_{k=0}^n c_k(r, s) P_k(i).$$ By the way, this so called Hamming scheme is an association scheme and their very useful property is (as you wrote) that the intersection numbers are the same for all pairs $(a,b)$ with $d(a,b)=k$. So, they only depend on the distance between $a$ and $b$ but not on the concrete codewords $a$ and $b$. There is a nice and short introduction of association schemes in "Handbook of Combinatorics Volume I". There you can also find how you can use these Krawtchuck polynomials and linear programming to compute an upper bound for the cardinality of error-correcting codes.
Proof of the Lemma: Without loss of generality, we can assume $x=(x_1,...,x_i,0,...,0)$ with $x_1,...,x_i\not=0$. Choose $k$ positions $p_1,...,p_k$ such that $$0<p_1<p_2<...<p_j\leq i < p_{j+1} <...<p_k\leq n.$$ For a fixed $j$, there are $\binom{i}{j}\binom{n-i}{k-j}$ possibilities to choose these positions. First, we choose the positions $p_1,...,p_j$ which are in the set $\{1,2,...,i\}$. There we have $\binom{i}{j}$ possibilities. Then, we choose the remaining $k-j$ positions $p_{j+1},...,p_{k}$ which are in $\{i,i+1,...,n\}$. Thus, $\binom{n-i}{k-j}$ possibilities.
Let $N$ be the set containing all codewords in $\mathbb{Z}_q^n$ of weight $k$ that have their nonzero entries in the positions $p_1,...,p_k$. Then, we can rewrite the left-hand side of (1) as follows $$\sum\limits_{y\in \mathbb{Z}_q^n\\wt(y)=k}\omega^{<x,y>}= \sum\limits_{j=0}^n \binom{i}{j}\binom{n-i}{n-j} \sum\limits_{y\in N}\omega^{<x,y>}.\qquad\qquad (2)$$ Now, we only have to look at the last sum which can be written as follows $$\sum\limits_{y\in N}\omega^{<x,y>} = \sum\limits_{y\in N}\omega^{x_1y_1+...+x_ny_n} =\sum\limits_{y_{p_1}\in\mathbb{Z}_q\setminus \{0\}}... \sum\limits_{y_{p_k}\in\mathbb{Z}_q\setminus \{0\}} \omega^{x_{p_1}y_{p_1}+...+x_{p_k}y_{p_k}}$$ since $y_l=0$ for all positions $l\not\in\{p_1,...,p_k\}$. We also know that $x_l=0$ for the $k-j$ positions $l\in \{p_{j+1},...p_k\}$. Therefore, we get $$\sum\limits_{y_{p_1}\in\mathbb{Z}_q\setminus \{0\}}... \sum\limits_{y_{p_k}\in\mathbb{Z}_q\setminus \{0\}} \omega^{x_{p_1}y_{p_1}+...+x_{p_k}y_{p_k}} = \left( \sum\limits_{y\in \mathbb{Z}_q\setminus\{0\}} \omega^{0\cdot y} \right)^{k-j} \prod\limits_{l=1}^j \sum\limits_{y\in \mathbb{Z}_q\setminus\{0\}} \omega^{x_{p_l}y}. $$ The first sum equals $|\mathbb{Z}_q\setminus\{0\}|=q-1$ and the second sum is $(-1)$ since $\omega^{x_{p_l} \cdot y}$ is a character of $\mathbb{Z}_q$. Thus, we get $$\left( \sum\limits_{y\in \mathbb{Z}_q\setminus\{0\}} \omega^{0\cdot y} \right)^{k-j} \prod\limits_{l=1}^j \sum\limits_{y\in \mathbb{Z}_q\setminus\{0\}} \omega^{x_{p_l}y} = (q-1)^{k-j} (-1)^j. $$ The lemma follows now by combining this with (2).