I assume that the encoding process for a $r$-th order Reed-Muller code
is the following.
Denote the $\sum_{i=0}^r\binom{m}{i}$ data bits as
$$\begin{align*}
&d_0,\\
&d_1, \quad d_2, \quad \ldots, \quad d_m,\\
&d_{1,2}, d_{1,3}, \ldots, d_{1,m},\quad d_{2,3}, d_{2,4},\ldots,
d_{2,m}, \quad\ldots,\quad d_{m-1, m}\\
&\vdots \qquad \ddots\\
&d_{1,2,\ldots, r}, \quad\ldots, \quad d_{m-r+1, m-r+2, \ldots, m},
\end{align*}$$
where the $i$-th row above lists $\binom{m}{i}$ data bits.
Form the polynomial
$$\begin{align*}
d(x_1, x_2, \ldots, x_m) &= \qquad d_0 \\
&\quad \oplus\quad d_1 x_1 \oplus d_2x_2 \oplus \cdots \oplus d_mx_m\\
&\quad \oplus \quad d_{1,2}x_1x_2 \oplus d_{1.3}x_1x_3
\oplus \cdots d_{m-1,m}x_{m-1}x_m\\
&\quad \oplus \quad \cdots \\
&\quad \oplus \quad d_{1,2,\ldots,r} x_1x_2\cdots x_r
\oplus \cdots d_{m-r+1,m-r+2,\ldots,m}x_{m-r+1}x_{m-r+2}\cdots x_m
\end{align*}$$
of degree $r$ in $m$ binary variables $x_1, x_2, \ldots, x_m$. Note that
the $i$-th row has all the $\binom{m}{i}$ terms of degree $i$.
Transmit (or record) the $2^m$ values of this polynomial
(Boolean function) as the
codeword. That is,
$$\mathbf{d} = \left(d(0,0,\ldots,0), d(0,0,\ldots, 0,1),
d(0, 0, \ldots, 1, 0),\ldots, d(1,1, \ldots, 1)\right) $$
is transmitted. Note that this just the truth table for this Boolean function.
The Hamming weight of $(d(x_1, x_2, \ldots, x_m)$ is the
number of entries in $\mathbf d$ that have value $1$. In what
follows, we will need the following property: The Hamming weight
of any polynomial is odd if and only if the polynomial is of
degree $m$, that is, $d_{1,2,\ldots, m} = 1$ and thus the
term $x_1x_2\cdots x_m$ is included in the polynomial
Given the codeword (or truth table) $\mathbf d$, how can we determine the
value of the data bit $d_{1,2,\ldots,r}$ which is the
coefficient of the term $x_1x_2\cdots x_r$ in $d(x_1,x_2,\ldots, x_m)$?
Well, let us consider the polynomial
$$g(x_1,x_2,\ldots, x_m)=x_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)$$
where we assume that for all $i, 1 \leq i \leq m$, all occurrences of $x_i^2$ on the
right side of the above equation have been replaced
by $x_i$ since $x_i$ takes on values $0$ and $1$ only. Since
$d(x_1, \ldots, x_m)$ is of degree $r$ of the form
shown above, we conclude that
$x_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)$ includes the
term $x_1x_2\cdots x_m$ if and only if $d_{1,2,\ldots,r} = 1$,
that is, $d(x_1, \ldots, x_m)$ contains the term $x_1x_2\cdots x_r$.
It follows that the truth table of $x_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)$
contains an odd number of $1$ entries if and only if $d_{1,2,\ldots,r} = 1$.
Now, the truth table of $x_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)$
has $0$s wherever at least one of $x_{r+1}, x_{r+2}, \ldots, x_m$
has value $0$, and the nonzero entries, if any, must be in
the $2^{r}$ positions where
$x_{r+1}=x_{r+2} = \cdots = x_m = 1$ are fixed, and only
$x_1, x_2, \ldots, x_r$ are varying. Furthermore, the entries
in these $2^r$ positions
must be the same as the entries in $\mathbf d$ -- which we know.
Are there an odd number of $1$ entries in $\mathbf d$ in these $2^r$ positions?
We can find out by computing the parity (also known as the
Exclusive-OR sum) of these $2^r$ bits in $\mathbf d$. Formally, we can
write
$$\begin{align*}
d_{1,2,\ldots, r} = \quad &d(0,0,\ldots, 0, 0, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 0, 1, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 1, 0, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 1, 1, 1, 1, \ldots, 1)\\
\quad \oplus &\cdots\\
\quad \oplus &d(1,1,\ldots, 1, 1, 1, 1, \ldots, 1).\\
\end{align*}$$
Any reader who has not been totally bored out of his gourd
by the prolix presentation above will have realized that
similar parity sums can be set up for the coefficients of
all the other terms of degree $r$. All we have to do is
create a $g(x_1, x_2, \ldots, x_m)$ by multiplying
$d(x_1, x_2, \ldots, x_m)$ by all the missing variables
and then compute the parity of the appropriate $2^r$ bits
in $\mathbf d$. (To paint refined gold and gild the lily,
the missing variables were $x_{r+1}, x_{r+2},\ldots,x_m$ when
we were trying to determine the coefficient of $x_1x_2\cdots x_r$).
But, wait, there is more. Once we have determined the coefficients
of the degree $r$ terms, it is necessary to modify $\mathbf d$
to be the truth table not of the $d(x_1, x_2, \ldots, x_m)$ that
we began with but rather
of a polynomial $\hat{d}(x_1, x_2, \ldots, x_m)$ of degree $r-1$.
This new polynomial is just the terms of degree $r-1$ or less
in $d(x_1, x_2, \ldots, x_m)$, and the modification is just
subtracting off the truth table of the polynomial
$$d_{1,2,\ldots,r} x_1x_2\cdots x_r
\oplus \cdots d_{m-r+1,m-r+2,\ldots,m}x_{m-r+1}x_{m-r+2}\cdots x_m
$$
from $\mathbf d$ (same as Exclusive-ORing in the truth table into $\mathbf d$).
If this is not done, then multiplying by the missing variables
$x_r, x_{r+1}, \ldots x_m$ in attempting to calculate $d_{1,2,\ldots, r-1}$
will not work because the term $d_{1,2,\ldots,r} x_1x_2\cdots x_r$
will create additional $x_1x_2\cdots x_m$ terms that interfere
in the calculation of $d_{1,2,\ldots, r-1}$.
Thus, the decoding
works in stages: first determine the coefficients of the
degree-$r$ terms and subtract
them off, then determine the coefficients of the
degree-$(r-1)$ terms and subtract them off,
and continue until all the terms have been calculated.
But if you order today, I will throw in error-correction for
free! All of the above works if the transmitted (or recorded)
codeword $\mathbf d$ is received (or read out) without any
errors. But when some of the bits in $\mathbf d$ are incorrect,
we cannot be sure that the parity calculation
$$\begin{align*}
d_{1,2,\ldots, r} = \quad &d(0,0,\ldots, 0, 0, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 0, 1, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 1, 0, 1, 1, \ldots, 1)\\
\quad \oplus &d(0,0,\ldots, 1, 1, 1, 1, \ldots, 1)\\
\quad \oplus &\cdots\\
\quad \oplus &d(1,1,\ldots, 1, 1, 1, 1, \ldots, 1).\\
\end{align*}$$
will give the correct value of $d_{1,2,\ldots, r}$. After
all, an odd number of bits participating in the parity
calculation might be wrong. Thus, the above parity calculation
should be treated as an estimate of the value of $d_{1,2,\ldots, r}$.
Fortunately, we can obtain $2^{m-r} - 1$ other estimates from
of $d_{1,2,\ldots, r}$ from other bits in $\mathbf d$.
For example, the polynomial
$$\begin{align*}
g(x_1,x_2,\ldots, x_m)
&= \bar{x}_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)\\
&= (x_{r+1} \oplus 1)x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)\\
&= x_{r+1}x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)
\oplus x_{r+2}\cdots x_md(x_1,x_2,\ldots, x_m)
\end{align*}$$
has degree $m$ if and only if $d_{1,2,\ldots, r} = 1$
but its truth table has nonzero entries in $2^r$ other
locations, viz. where $x_{r+1}=0, x_{r+2}= x_{r+3} = \cdots = x_m = 1$.
In this fashion, we do $2^{m-r}$ different parity calculations
for $d_{1,2,\ldots, r}$ using all the $2^m$ bits in $\mathbf d$,
and then take a majority vote of the different estimates
to estimate the value of $d_{1,2,\ldots, r}$. As long as no
more that $2^{m-r-1}-1$ bits in the received $\mathbf d$
are incorrect, the majority vote results in a correct estimate,
that is, up to $2^{m-r-1}-1$ bit errors can be corrected:
a tie vote indicates that a detectable but uncorrectable
error pattern has occurred.
The rest of the decoding proceeds as described before with the
obvious changes of doing multiple parity calculations in estimating
each data bit, and what is subtracted off between stages is
the estimated terms of that degree. Note that more parity
calculations can be made in later stages, which can make
assurance doubly sure if the previous stages resulted in
correct decoding, but provide no additional error protection
if incorrect decoding occurred at an earlier stage.
In summary, Reed-Muller decoding is a lot like peeling
an onion. It has to be done one layer at a time, and inexperienced
decoders (and cooks) tend to weep a lot.