$$
\newcommand{\bra}[1]{\left<#1\right|}\newcommand{\ket}[1]{\left|#1\right>}\newcommand{\bk}[2]{\left<#1\middle|#2\right>}\newcommand{\bke}[3]{\left<#1\middle|#2\middle|#3\right>}\newcommand{\proj}[1]{\left|#1\right\rangle\left\langle#1\right|}
$$
Much of the functionality here is the same as the Bernstien-Vazirani algorithm, if that helps. The following is more or less copy and pasted from some lecture notes I prepared at some point. It explains it in a slightly different way to the direction you're coming at it from, but hopefully gets you going in the right direction.
The circuit is, in principle, the same as for the Bernstein-Vazirani Algorithm, except that since the output of the function evaluation is $n$ bits, the second register, which is used for the reversible function evaluation, is also $n$ bits,
This function evaluation is represented by the controlled-controlled...-controlled-U gate, acting as $\ket{x}\ket{y}\mapsto\ket{x}\ket{y\oplus f(x)}$.
We start with the first set of Hadamards producing an equally weighted superposition of all strings of $n$ bits,
$$
\ket{0}\ket{0}\rightarrow\frac{1}{\sqrt{2^n}}\sum_{x\in\{0,1\}^n}\ket{x}\ket{0}.
$$
This is followed by the function evaluation,
$$
\rightarrow\frac{1}{\sqrt{2^n}}\sum_{x\in\{0,1\}^n}\ket{x}\ket{f(x)}.
$$
At this point, you could measure the second register. This will return a random value of $f(r)$, so that the overall state of the system is in
$$
\frac{1}{\sqrt{2}}(\ket{r}+\ket{r\oplus s})\ket{f(r)}.
$$
So, the point is that you don't have to actively prepare different values of $r$; they will be selected for you by the measurement. We'll have to make a bit of an argument later on about how likely it is that we get new information each time we make such a measurement, but it'll all work out.
Actually, I don't usually think about performing the measurement at this point; it's unnecessary (but would allow you to avoid the density matrix formalism in the following calculation). Instead, calculate the action of the final set of Hadamards, yielding the final output state
$$
\frac{1}{2^n}\sum_{x,z}(-1)^{x\cdot z}\ket{z}\ket{f(x)}.
$$
Here,
$$
x\cdot z= x_1z_1\oplus x_2z_2\oplus x_3z_3 \oplus\ldots \oplus x_nz_n,
$$
where $x_k$ is the $k^{th}$ bit of $x$.
Now we collect unique values of $f(x)$,
$$
\frac{1}{2^n}\sum_{z,f(x)}\left((-1)^{x\cdot z}+(-1)^{(x\oplus s)\cdot z}\right)\ket{z}\ket{f(x)}.
$$
One can therefore verify that the output of the algorithm is
$$
\frac{1}{2^{n-1}}\sum_{z: s\cdot z=0}(-1)^{x\cdot z}\ket{z}\sum_{f(x)}\ket{f(x)}.
$$
The state of the first register clearly contains the information about $s$, but we need to extract it. Measuring the second register is not necessary (we could do it, but it doesn't help). Instead, let's trace out the second register, so we get
$$
\rho=\frac{1}{2^{2n-2}}\sum_{z: s\cdot z=0}\sum_{y: s\cdot y=0}\sum_x(-1)^{x\cdot(y\oplus z)}\ket{y}\bra{z}.
$$
By performing the sum over $x$, we are left with
$$
\rho=\frac{1}{2^{n-1}}\sum_{s\cdot z=0}\proj{z}.
$$
Measurement of the first register yields a binary string $z$ where $s\cdot z=0$. If we had $n-1$ such examples which are linearly independent, we would be able to determine $s$. This requires repeated application of the algorithm to find enough strings (this is the 'Fourier Sampling' part).
We must now justify that a linear number of applications is sufficient to find enough strings with high probability. In the absolute worst case, when we have found $n-2$ vectors, we must find the one remaining bit of information when there are still $2^{n-1}$ strings $s\cdot z=0$ to sample. These must constitute $1-2^{n-2}/2^{n-1}=\frac12$ of the space. Hence, the average number of trials to find one of these vectors is given by
$$
\frac{1}{2}\sum_{n=0}^\infty\frac{n+1}{2^n},
$$
which is readily evaluated using the following identities,
\begin{eqnarray}
\sum_{n=0}^\infty r^n&=&\frac{1}{1-r} \nonumber\\
\frac{d}{dr}r\sum_{n=0}^\infty r^n&=&\sum_{n=0}^\infty (n+1)r^n. \nonumber
\end{eqnarray}
On average, we always find another linearly independent vector by making $2$ samples, i.e. within $2n$ steps, we have a high probability of determining $s$. It is this final part of the argument, requiring some classical post-processing that differentiates Simon's algorithm from the Bernstein-Vazirani algorithm. Essentially, the difference comes from not knowing the eigenvectors of $U$ in advance. If we did, we could prepare the second register in a fixed state. Instead, it is prepared in a superposition of different eigenstates (which coincides with a nice state to prepare) and we have to rely on a certain amount of randomness to sample the elements we need.
