Why are the $p$-oldforms $f(z)$ and $f(pz)$ linearly independent at level $\Gamma_0(pN)$? Let $f$ be a newform (normalized eigenform) of weight $k$ and level $\Gamma_0(N)$. Fix $p$ not dividing $N$ and set $f_p(z)=f(pz)$. 

Viewing $f$ and $f_p$ at level $\Gamma_0(pN)$, why are they linearly independent? 

Here is what I'm thinking: Suppose there is some nontrivial relation 
$$
af+bf_p=0.
\tag{1}$$ The Hecke eigenvalues for $f$ are the same at level $pN$, except possibly at $p$. Thus, for any $\ell\neq p$, applying $T_\ell-a_\ell(f)$ to equation (1) shows that $T_\ell f_p=a_\ell(f)f_p$. Thus, $f$ and $f_p$ have the same Hecke eigenvalues outside $p$. 
Now, at this point, I have heard mumblings that the "multiplicity one theorem" shows that $f=f_p$, but I can't seem to find a suitable reference. Simply Googling "multiplicity one" in this context leads to many high-brow results... Diamond and Shurman do not state the multiplicity one theorm in their book (my main reference thus far), but they do cite Miyake's book as a reference. I looked up the theorem there (I believe it is Theorem 4.6.12), and it seems to require that $f$ and $f_p$ both be of level $N$, which (I am fairly certain) they are not. Can someone clarify this for me and possibly provide a good reference for the classical multiplicity one theorem? 

From the Fourier series it is obvious that $f \ne cf(pz)$. 
$$f(z)=\sum_n a(n) e^{2i\pi n z},\qquad f(pz)=\sum_n a(n) e^{2i\pi np z}=\sum_n a(\frac{n}p)1_{p|n} e^{2i\pi n z}$$
The sequence $a(n)$ is not a scalar multiple of $a(\frac{n}p)1_{p|n}$.
The Fourier series is uniquely defined by the function : $$a(n)=\int_{iy}^{iy+1}f(z)e^{-2i\pi nz}dz$$
