This question already has an answer here:
I will prove your edited formula, which is much nicer.
There is a simple solution (two lines) via generating functions. In particular, read the 2nd example which uses convolutions.
Consider the generating functions which gives us $ a_n = { 2n \choose n } $. Let this be denoted by $ a(x) $.
Consider the generating functions which gives us $ b_ n = 4 ^n $. Clearly, this is generated by $ b(x) = \frac{1}{1-4x} $.
Then, the statement claims that
$$ a(x) ^ 2 = b( x). $$
All that is left to do, is to show that
$$ a(x) = \frac{1}{ \sqrt{ 1 – 4 x } } $$
This is obvious / well known.
Here’s an elementary combinatorial proof. Let
$$q(m)=\sum_{k=0}^m{2k\choose k}{2m-2k\choose m-k}.$$
As in the question, $q(m)$ counts the number of times a monotonic $m$-path (by which we mean a lattice path from $(0,0)$ to $(m,m)$ consisting of unit steps to the east and unit steps to the north) touches the diagonal (including at its beginning and its end), summed over all paths. We will prove that $q(m)=4^m$ for all $m$. It clearly is for $m=0$. Suppose $q(m)=4^m$ is true for all $m$ up to $n-1$, where $n\geqslant1$.
The induction hypothesis allows us to count the number of times a monotonic $n$-path touches the diagonal in a cleaner way: Any path touches the diagonal in the origin. For a path $p$ denote $t_p$ the x (or y) coordinate of the next meeting point with the diagonal. (So $t_p\in\{1,\ldots,n\}$ for any $n$-path $p$.) The number of paths from $(0,0)$ to $(k,k)$ with $t_p=k$ is $2C_{k-1}$ (see Wikipedia; note that in our case the path should not only stay at the same side of the diagonal but it shouldn’t even touch it, which is why we have $C_{k-1}$ and not just $C_k$). Here as usual $C_n=\frac1{n+1}{2n\choose n}$ denotes the $n$th Catalan number.
A monotonic $n$-path from $(0,0)$ to $(n,n)$ with $t_p=k$ touches the diagonal once at the origin and then follows a monotonic $n-k$-path. For fixed $k$, the number of touches with the diagonal, counting only those $n$-paths that satisfy $t_p=k$, is $2C_{k-1}\cdot q(n-k)$ plus the number of possible paths with $t_p=k$ (because of that one touch in the origin), which is $2C_{k-1}\cdot\binom{2n-2k}{n-k}$. Summing over $k$ gives
$$q(n)=\sum_{k=1}^n2C_{k-1}\left(q(n-k)+\binom{2n-2k}{n-k}\right).$$
By the induction hypothesis,
$$q(n)=\sum_{k=1}^n2C_{k-1}\left(4^{n-k}+\color{red}{\binom{2n-2k}{n-k}}\right).$$
What follows is an algebraic manipulation to show that this is indeed equal to $4^n$.
Using the identity $2C_k=\frac 2{k+1}\binom{2k}{k}=\color{green}{4\binom{2k}{k}}-\color{blue}{\binom{2k+2}{k+1}}$ this can be rewritten as
$$\underbrace{\color{green}4\sum_{k=1}^n\color{green}{\binom{2k-2}{k-1}}\color{red}{\binom{2n-2k}{n-k}}}-\underbrace{\sum_{k=1}^n\color{blue}{\binom{2k}{k}}\color{red}{\binom{2n-2k}{n-k}}}+\underbrace{\sum_{k=1}^n\left(\color{green}{4\binom{2k-2}{k-1}}-\color{blue}{\binom{2k}{k}}\right)4^{n-k}}$$
which is
$$\underset{\text{(Ind. Hyp.}\\\text{for }m=n-1)}{\underbrace{4\cdot4^{n-1}}}-\underset{\text{(Def. of }q(n))}{\underbrace{q(n)+{2n\choose n}}}+\underset{\text{(Telescoping)}}{\underbrace{4^n-{2n\choose n}}}=2\cdot4^n-q(n).$$
Hence
$$q(n)=2\cdot4^n-q(n)$$
and the conclusion follows.
Consider the taylor series expansion
$$(1/(1-x))^\frac12=\sum_{k \ge 0} \frac{(2k-1)!!}{(2k)!!}x^k,\tag{1}$$
whose square is of course the desired sequence of all coefficients $1$ in front of the $x$ powers, and whose coefficients are successively $1,1/2,3/8,5/16,35/128,$ which match your values for $f(p^n).$ [I think it should be easy to manipulate the double factorial term format of the coefficients in $(1)$ to your version of $2n-1$ choose $n$ over $2^(2n-1).$] In $(1)$ the double factorials are as usual by successively decreasing terms by $2$, e.g. $5!!=5\cdot 3 \cdot 1.$ The explicit form of $(1)$ in terms of the double factorial ratios for coefficients comes from a formula in Gradshteyn & Ryzhic, which in the 1980 edition of “Tables of Integrals, Series, and Products” appears on page 21, and is formula 1.112 (4) with its $x$ replaced by $-x.$
Anyway given that the expansion $(1)$ is the correct one for that squareroot, and also given that its coefficients match up with the ones you have conjectured, it follows by the way ordinary series multiply that the terms are correct for the Dirichlet square of that part of the product for powers of $p$ to wind up as the ordinary series of $1$ in front of each power.
I realize I’m missing some details as to setting this all up. But I think it interesting that the usual taylor series of the squareroot of $1/(1-x)$ has its coefficients matching your $f(p^n)$, and in one book on analytic number theory I have consulted, such switching back and forth between Dirichlet series and ordinary ones is justified when the given function is prime-independent and multiplicative. I’ll add more on this eventually but wanted to enter the answer now in case it’s of any immediate use.