To simplify the notation of the formulas in the frequency domain we $z = exp(j \omega T_s)$ where $T_s$ is the sampling period, and $\omega$ some angular frequency, under these conditions that discrete Fourier transform coincides with z transform .
If you assume $F(z)$ to have a short input response and to change slowly enough that you can use it as a time invariant filter in a block of $N$ samples.
You can compute the energy in the frequency domain as $m(z) \bar{m}(z)$$m(z) m(z)^*$, I am using the $*$ suffix to denote complex conjugation.
$$ \begin{eqnarray} m(z)^*m(z) &=& (F(z) k(z) + s(z))^*(F(z) k(z) + s(z)) \\ &=& \underbrace{|F(z)k(z)|^2}_{\textrm{echo energy}} + 2 Re(F(z)k(z)s(z)*) + \underbrace{|s(z)|}_{\textrm{external sound}} \end{eqnarray} $$$$ \begin{eqnarray} m(z)^*m(z) &=& (F(z) k(z) + s(z))^*(F(z) k(z) + s(z)) \\ &=& \underbrace{|F(z)k(z)|^2}_{\textrm{echo energy}} + 2 Re(F(z)k(z)s(z)^*) + \underbrace{|s(z)|}_{\textrm{external sound}} \end{eqnarray} $$
And you can increase alpha for a more stable estimate, or decrease it to have a faster estimate.
Notes on complex arithmetic
To answer the question about the term 2 Re(.) left on the comment. The short answer is this is the crossed terms in the expansion of the product in the previous.
I skipped some steps there to avoid making long arrays of equations since this is a very common result in complex manipulations.
$$\begin{eqnarray} (x + y)^* (x + y) &=& (x^* + y^*)(x + y) \\ &=& (x^* x + y^* y + x^*y + y^*x \\ &=& x^* x + y^*y + (x^*y) + (x^*y)^* \end{eqnarray}$$
Let $x = a+bj$, with $j^2=-1$, term $$x^*x=|x|^2$,
$$\begin{eqnarray}x^*x &=& (a-bj)(a+bj) \\&=* (a^2 - b^2j^2 -baj +abj) \\&=& (a^2 -b^2(-1)) \\&=& (a^2 + b^2) \\&=& |x|^2\end{eqnarray}$$
Similarly, $y^*y = |y|^2$.
The terms $x^*y + y^*x$ are conjugate from each other, thus they could be written as $w^* + w$, since $w = Re(w) + j Im(w)$ and $w^* = Re(w) - j Im(w)$, $w + w^* = 2 Re(w)$, where $Re(w)$ denotes the real part, and $Im(w)$ denotes te imaginary part.
Why the energy is the square of the amplitude
I will explain this using a Mass-spring system, let a mass $m$ attached to a sprint of constant $k$ oscillating with amplitude $A$. The maximum potential energy of the spring is achieved when it is deformed by $A$, and can be calculated as $\int_{0}^{A} k \, x\, dx = k A^2 / 2$. Also, if there are multiple oscilation modes (different frequencies), each one conserves energy separately, this is why we can talk about power at specific frequencies.