Non standard solution to $f(x) = \frac{1}{2}\Big(f(\frac{x}{2}) + f(\frac{1+x}{2})\Big)$ Final update: on 11/29/2019: I have worked on this a bit more, and wrote an article summarizing all the main findings. You can read it here.
This functional equation appears in the following context. Let $\alpha\in[0,1]$ be an irrational number (called seed) and consider the sequence $x_n=\{2^n \alpha\}$. Here the brackets represent the fractional part function. In particular, $\lfloor 2x_n\rfloor$ is the $n$-th digit of $\alpha$ in base $2$.
The values $x_n$ are distributed in a certain way due to the ergodicity of the underlying process. The density associated with this distribution is the function $f$, and for the immense majority of seeds $\alpha$ that density is uniform on $[0, 1]$, that is, $f(x) = 1, x \in [0, 1]$. Such seeds $\alpha$ producing the uniform density are sometimes called normal numbers; their digit distribution is also uniform.  
However, the functional equation $f(x) = \frac{1}{2}\Big(f(\frac{x}{2}) + f(\frac{1+x}{2})\Big)$ may have plenty of other solutions. Such solutions are called non-standard solutions. Can you find a seed $\alpha$ producing a non-standard solution, with an explicit form for $f$? Maybe a step-wise uniform function? The set of seeds producing non-standard solutions is known to have Lebesgue measure zero, but there are infinitely many such seeds. 
All rational seeds $\alpha$ work, but they produce a discrete distribution. Thus their density is of the discrete type. Here however, I am interested in a continuous function $f$, even if it has infinitely many points of discontinuity (that is, a function $f$ continuous almost everywhere: the set of discontinuity points has Lebesgue measure zero.) 
Update
I am looking for a function $f$ that is a density on $[0, 1]$, so there are additional constraints here: $\int_0^1 f(x)\,dx = 1$ and $f(x) \geq 0$. However, note that if $f$ is a solution, then $cf$ is also a solution regardless of the constant $c$. So any solution can be normalized to integrate to one. Also, $cf+d$ is also a solution ($c, d$ constants).
Second update
Below is a density satisfying all the requirements. Actually, the plot below represents its percentile distribution. It was produced with a seed $\alpha$ built as follows: its $n$-th binary digit is $1$ if $\mbox{Rand}(n) < 0.75$, and $0$ otherwise, using a pseudo random number generator. Note that $P._{25} = 0.5$ and corresponds to a dip ($P._{25}$ denotes the $25$-th percentile.) Dips are everywhere, only the big ones are visible. By contrast, the percentile distribution for the uniform case (if you replace $0.75$ by $0.50$ in $\mbox{Rand}(n) < 0.75$) is a straight line, with no dips. 

Note: I eventually answered my question, see the second answer.
 A: Not an answer, just some thoughts.
Let's try applying Fourier transform:
$$f(x)=\frac{1}{\sqrt{2 \pi}}\int_{-\infty}^\infty e^{i k x} g(k) dk$$
The equation becomes:
$$\int_{-\infty}^\infty \left(2 e^{i k x}-e^{i k x/2} (1+e^{i k /2}) \right) g(k) dk=0$$
Making a simple substitution:
$$\int_{-\infty}^\infty e^{i k x} \left(g(k)- (1+e^{i k }) g(2k) \right)  dk=0$$
We can consider a particular case (or is this the general case?):
$$g(k)=(1+e^{i k }) g(2k)$$
Examples:
$$g(0)=0$$
For any $a \neq 0$:
$$g(a)=C$$
$$g(2a)= \frac{C}{1+e^{ia}}$$
$$g(2^n a)= C \prod_{q=1}^{n-1} \frac{1}{1+e^{i q a}}$$
In the same way we can find $g(2^{-n} a)$.
For $a$ not being a rational multiple of $\pi$ the infinite product doesn't converge, and the function itself is not nice.
If we call:
$$G_n= -\log g(2^n), \qquad g(1)=1$$
Then the plot $G_n$ looks like this:

Not sure if this could lead to a non-trivial solution to the original equation though.
A: Here I discuss in more detail the case mentioned in the section "second update". This special distribution was produced with a seed $\alpha$ built as follows: its $n$-th binary digit is $1$ if $\mbox{Rand}(n)< p$, and $0$ otherwise, using a pseudo random number generator. I used $p=0.75$ in my example. 
Now $x_n$ (introduced in the first paragraph in my question) is a random variable, and we have:
$$x_n=\sum_{k=1}^\infty \frac{d_{n+k}}{b^k}.$$
Here $b$ is the base ($b=2$), and $d_{n+k}$ is the $(n+k)$-th digit of $\alpha$ in base $b$. Furthermore, by construction, these digits are identically and independently distributed with a Bernouilli distribution of parameter $p$. Thus, using the convolution theorem, the characteristic function of $x_n$ is
$$\phi(t; p, b) = \prod_{k=1}^\infty \Big(1-p(1-\exp \frac{it}{b^k})\Big).$$
Take the derivative of the inverse Fourier transform (see section inverse formula here) and you obtain
$$f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-ity} \phi(t; p,b) dy.$$
If $p=0.5$ and $b=2$ we are back to the uniform case. If $p\neq 0.5$ then the solution is quite special: the density $f$ is nowhere continuous it seems. See picture below for $p=0.55, b=2$.

Now we should prove that this case is ergodic, for the functional equation to apply. I also tried to check with some sampled values of $x$ to see whether  $f(x) = \frac{1}{2}\Big(f(\frac{x}{2}) + f(\frac{1+x}{2})\Big)$, but the function being discontinuous everywhere, and since I got its value approximated probably to no more than two decimals, it is not easy. The distribution attached to this density has the following moments: 


*

*Expectation: $\frac{p}{b-1}.$ 

*Variance: $\frac{p(1-p)}{b^2-1}.$
An Excel spreadsheet with the computation of $x_n$ up to $n=200,000$ and with a precision of 14 decimals, is available upon request. You can interactively change $b$ or $p$ and see the result on the chart. Despite dealing with more than 200K digits, the computations are done very efficiently. Finally, the functional equation can be adapted to base $b$, provided $b$ is an integer. It becomes:
$$f(x) = \frac{1}{b}\sum_{m=0}^{b-1}f\Big(\frac{x+m}{b}\Big).$$
A: We are told that $f\in L^1([0,1])$ and that the equation $$
f(x)=\frac 1 2 \left( f\left(\frac x 2 \right)+f\left(\frac{1+x} 2\right)\right)\tag{*}$$ holds, for almost all $t\in[0,1]$.
Let $a_n=\int_0^1\exp(2\pi i n x) f(x)dx$.   Formula (*) gives:
$$ a_n = \frac1 2 \int_0^1 \exp(2\pi i n x) f(x/2)dx + \frac 1 2 \int_0^1 \exp(2\pi i n x)f((1+x)/2)  dx$$
$$ = \int_0^{1/2}\exp(2\pi i 2u ) f(u) du + \int_{1/2}^1 \exp(2\pi i 2u ) f(u) du $$ by substituting $u=2x$ and so
$$a_n = \int_0^1 \exp(2\pi i n 2u) f(u)du = a_{2n},$$
and so $a_n=a_{2n}=a_{4n}=a_{8n}\cdots.$
If $n\ne0$ we have $a_{2^kn} = a_n$ for all $k$, so for such an $n$,
$$a_n = \lim_{k\to\infty} a_{2^kn} = 0$$
by the Riemann Lebesgue lemma. This implies $f$ is a constant, (the constant $a_0$) almost everywhere.
