I looked at that paper of Fejér, and although he does discuss some related things, I did not find that exact statement in his paper (granted, my German is probably not much better than yours). So, I will just provide a proof with the only hint from the Fejér's paper being that we have to use Cauchy--Schwarz inequality somewhere.
In fact, I will show that $r_N = \left(1 - \frac{1}{N}\right)$ works for all such functions, so it does not even depend on $f$ (proof also works for $r_N = \left(1 - \frac{c}{N}\right)$ for any fixed positive $c$).
So, we have to bound
$$\sum_{n = 0}^N (1 - r_N^n) a_n e^{i\theta n} + \sum_{n = N+1}^\infty r_N^n a_n e^{i\theta n}.$$
Fix some big number $K$ and split the first sum additionally into the terms $n \le K$ and $n > K$, where $K$ is chosen so that $\sum_{n > K} n |a_n|^2 < \varepsilon$. This way we will have $I_1 + I_2 + I_3$. We will show that for big enough $N$ each of these terms is at most $\sqrt{\varepsilon}$, this will be enough since $\varepsilon > 0$ can be arbitrary.
For $I_1$ there is nothing to prove, as long as $K$ is fixed it tends to $0$ uniformly since it has $K$ terms each of which tends to $0$ in absolute value regardless of $\theta$. For $I_2$ and $I_3$ we will use the Cauchy--Schwarz inequality. We start with $I_2$:
$$|I_2| \le \sum_{n = K+1}^N (1 - r_N^n) |a_n| \le \sum_{n=K+1}^N \frac{n}{N} |a_n| \le \sqrt{\varepsilon} \sqrt{\sum_{n = K+1}^N \frac{n}{N^2}} \le \sqrt{\varepsilon},$$
where in the first step we put absolute value into the sum, in the second step we used Bernoulli's inequality $(1+s)^k \ge 1 + ks$ for $s \ge -1, k \in \mathbb{N}$, in the third step we used Cauchy--Schwarz inequality, and in the last step we just used that each term in the sum is at most $\frac{1}{N}$ and the number of terms in the sum is at most $N$.
Now, we turn to $I_3$. Note that $r_N^{N} = \left (1 - \frac{1}{N}\right)^{N} \le \frac{1}{2}$ for big enough $N$ (in fact, this sequence tends to $\frac{1}{e}$). We will now use one more decomposition for $I_3$, denoting by $J_k$ the sum from $kN + 1$ to $(k+1)N$. On each such interval $r_N^n$ is at most $2^{-k}$. So, we have
$$|J_k| \le 2^{-k}\sum_{n = kN+1}^{(k+1)N} |a_n| \le 2^{-k}\sqrt{\varepsilon} \sqrt{\sum_{n = kN+1}^{(k+1)N} \frac{1}{n}} \le 2^{-k}\sqrt{\varepsilon},$$
where again the first step is putting absolute values, the second step is Cauchy--Scwarz inequality, and in the last step there are $N$ terms in the sum, each of which is at most $\frac{1}{N}$.
Summing this over $k$ we get that $|I_3| \le \sqrt{\varepsilon}$.