It is genuinely the case that $f(n) = \textrm{A007306}(n+1)$.
To start with, let's take the observation
Except for the initial one, this is the odd bisection of A002487. - Franklin T. Adams-Watters, May 25 2015
We're quite happy to discard the initial one, because that solves the offset problem! A002487 is Stern's diatomic sequence and has the simple recurrence $$\begin{eqnarray}\textrm{A002487}(0) & = & 0 \\
\textrm{A002487}(1) & = & 1 \\
\textrm{A002487}(2n) & = & \textrm{A002487}(n) \\
\textrm{A002487}(2n+1) & = & \textrm{A002487}(n) + \textrm{A002487}(n+1)\end{eqnarray}$$
Since we want to bisect it, let's define $N(k) = \textrm{A002487}(2k)$ and $D(k) = \textrm{A002487}(2k+1)$ (for numerator and denominator respectively). Then we have $$\begin{eqnarray}N(0) & = & 0\\
D(0) & = & 1\\
N(2k) & = & N(k) \\
D(2k) & = & N(k) + D(k) \\
N(2k+1) & = & D(k) \\
D(2k+1) & = & D(k) + N(k+1)
\end{eqnarray}$$
Borrowing another observation from OEIS,
$(\textrm{A002487}(n-1) + \textrm{A002487}(n+1))/\textrm{A002487}(n) = \textrm{A037227}(n)$ for $n \ge 1$. - Peter Bala, Feb 07 2017
and since $\textrm{A037227}(2k+1) = 1$ for all $k$, we have $N(k) + N(k+1) = D(k)$ and so we can rewrite that last case as $$D(2k+1) = 2D(k) - N(k)$$
Now, let's consider $f(n)$ as defined in the OP. Let $S_n$ be the set of distinct values obtained in step 4. The subsequences which generate $S_{2k+i}$ with $i \in \{0,1\}$ are the subsequences which generate $S_k$ and each of those subsequences with an $i$ appended, and so we have an iterative definition: $$\begin{eqnarray}S_0 & = & \{0\} \\
S_{2k} & = & S_k \cup \{2x \mid x \in S_k\} \\
S_{2k+1} & = & S_k \cup \{2x+1 \mid x \in S_k\}
\end{eqnarray}$$
We can partition $S_n$ into four subsets and hence partition $f(n) = |S_n|$ into a four-tuple $(a_n, b_n, c_n, d_n)$ where $$\begin{eqnarray}
a_n = |\{x \mid x \in S_n \wedge 2x \in S_n \wedge 2x+1 \not\in S_n \}| \\
b_n = |\{x \mid x \in S_n \wedge 2x \not\in S_n \wedge 2x+1 \in S_n \}| \\
c_n = |\{x \mid x \in S_n \wedge 2x \in S_n \wedge 2x+1 \in S_n \}| \\
d_n = |\{x \mid x \in S_n \wedge 2x \not\in S_n \wedge 2x+1 \not\in S_n \}| \\
\end{eqnarray}$$
Then we get a simple recurrence $$\begin{eqnarray}(a_0, b_0, c_0, d_0) & = & (1, 0, 0, 0) \\
(a_{2k}, b_{2k}, c_{2k}, d_{2k}) & = & (a_k + d_k, 0, b_k + c_k, b_k + d_k) \\
(a_{2k+1}, b_{2k+1}, c_{2k+1}, d_{2k+1}) & = & (0, b_k + d_k, a_k + c_k, a_k + d_k)
\end{eqnarray}$$
Clearly $c_k = d_k$ for all $k$, which simplifies things slightly.
Define $\nu_k = b_k + c_k$ and $\delta_k = a_k + b_k + 2 c_k = f(k)$. Then we have $$\begin{eqnarray}\nu_0 &= & 0 \\
\delta_0 & = & 1 \\
\nu_{2k} & = & b_k + c_k = \nu_k \\
\delta_{2k} & = & a_k + 2b_k + 3c_k = \nu_k + \delta_k \\
\nu_{2k+1} & = & a_k + b_k + 2c_k = \delta_k \\
\delta_{2k+1} & = & 2a_k + b_k + 3c_k = 2 \delta_k - \nu_k
\end{eqnarray}$$
Thus $(\nu_k, \delta_k)$ follows the same recurrence as $N(k), D(k)$ from the same initial conditions, and hence is identical. QED.
A more combinatorial proof would be nice. I wonder whether an alternative approach can be found in Calkin and Wilf's alternate bit set theorem, which certainly seems conceptually linked.
