If I have a large number of identical systems in identical quantum state $\Psi$ and an observable $A$ whose eigenstates are $\alpha_n$:
$$ \Psi = \sum_n c_n \alpha_n $$
I can get absolute values of $|c_n|^2$ by measuring $A$, but with which experiment I can get relative phases of $c_n$?
The phase ambiguity is a bit worse than you think. There is a global phase ambiguity in $|\Psi⟩$, for sure, but if the state $$ |\Psi⟩=\sum_n c_n |\alpha_n⟩ \tag 1 $$ is all you have around, then there is also a phase ambiguity in the phase of each individual $c_n$. This is because if you change $|\alpha_n⟩$ to $|\alpha'_n⟩=e^{i\theta_n}|\alpha_n⟩$, the transformed state is also an eigenstate of $A$ with eigenstate $|\alpha_n⟩$, so the two have exactly equal standing. If your universe only ever contains $|\Psi⟩$ (and a measurement device for $A$), then you only ever care about the absolute values $|c_n|^2$ of the coefficients.
The phase of these coefficients comes into play if you have two or more states, i.e. if you introduce some second state $$ |\Phi⟩=\sum_n b_n |\alpha_n⟩ \tag 2 $$ into the mix. Here transforming the eigenstates by some phase, as before, will still change the phases of the individual coefficients, but it will leave their differences $\arg(c_n)-\arg(b_n)$ unchanged. These are the phase differences that can be detected by experiment.
So, to come to the question, how do you detect them? There are a bunch of ways, but they all involve interference in some way. As a simple example, take your two states $|\Phi⟩$ and $|\Psi⟩$, and a superposition of them, $$ \frac{|\Phi⟩+|\Psi⟩}{\sqrt 2} = \sum_n \frac{c_n+b_n}{\sqrt 2} |\alpha_n⟩,\tag 3 $$ and measure $A$. Then you will measure the eigenvalue $\alpha_n$ with a probability $$P_n=\frac12 |c_n+b_n|^2$$ which is sensitive to the relative phase of the two coefficients.
Let's simplify your quantum state and talk about a qubit state:
$$| \psi \rangle = \cos\left(\frac{\theta}{2}\right) |{0}\rangle + e^{-i\phi}\sin\left(\frac{\theta}{2}\right) |{1}\rangle$$
In this case, the operator $A=Z$ the $\sigma^z$ Pauli matrix whose eigen-states are indeed $|0\rangle$ and $|1\rangle$. The relative phase $\phi$ between components is computationally important (Represent the state in the Bloch sphere). However, if you have only $Z$ operator in the world, you can only measure the angle $\theta$ without knowing any information about $\phi$. To determine $\phi$, it is necessary to rotate the state in the Bloch sphere or simply mix the amplitudes of the state. Namely, we need both $R_y\left(\frac{\pi}{2}\right)$ and $R_x\left(\frac{\pi}{2}\right)$: $$R_y\left(\frac{\pi}{2}\right)| \psi \rangle = \frac{1}{2} \left( \cos\left(\frac{\theta}{2}\right) + e^{-i\phi}\sin\left(\frac{\theta}{2}\right) \right) |{0}\rangle + \frac{1}{2} \left( \cos\left(\frac{\theta}{2}\right) - e^{-i\phi}\sin\left(\frac{\theta}{2}\right) \right) |{1}\rangle$$
In this case, the probability of obtaining $|0\rangle$ by measuring $Z$ is : $$p_0 = \frac{1}{2}+\cos(\phi)\cos\left(\frac{\theta}{2}\right)\sin\left(\frac{\theta}{2}\right)$$ and the probability of obtaining $|1\rangle$ by measurering $Z$ is : $$p_1 = \frac{1}{2}-\cos(\phi)\cos\left(\frac{\theta}{2}\right)\sin\left(\frac{\theta}{2}\right)$$
Determining $\cos(\phi)$ is not enough. Hence, it is necessary to apply $R_x\left(\frac{\pi}{2}\right)$:
$$R_x\left(\frac{\pi}{2}\right)| \psi \rangle = \frac{1}{2} \left( \cos\left(\frac{\theta}{2}\right) + ie^{-i\phi}\sin\left(\frac{\theta}{2}\right) \right) |{0}\rangle + \frac{1}{2} \left( \cos\left(\frac{\theta}{2}\right) - ie^{-i\phi}\sin\left(\frac{\theta}{2}\right) \right) |{1}\rangle$$.
In this case, the probability of obtaining $|0\rangle$ by measuring $Z$ is : $$p_0 = \frac{1}{2}+\sin(\phi)\cos\left(\frac{\theta}{2}\right)\sin\left(\frac{\theta}{2}\right)$$ and the probability of obtaining $|1\rangle$ by measurering $Z$ is : $$p_1 = \frac{1}{2}-\sin(\phi)\cos\left(\frac{\theta}{2}\right)\sin\left(\frac{\theta}{2}\right)$$
By this way, the relative phase angle $\phi$ is determined.
In general, the quantum state is multi-dimensional in Hilbert space: $$ |\Psi\rangle = \sum_{i=0}^{2^N - 1} \alpha_i |i\rangle$$
But, it is possible to write that state in the form of two ortho-normal states: For any eigenstate $|k\rangle$,
$$ |\Psi\rangle = \alpha_k |k\rangle + \bar{\alpha}_k |\bar{k}\rangle$$ such that: $$ |\bar{k}\rangle = \frac{1}{\sum_{i\neq k} |\alpha_i|^2} \sum_{i\neq k} \alpha_i |i\rangle$$ $$\bar{\alpha}_k = \sum_{i\neq k} |\alpha_i|^2$$
Now, we can redo the previous protocol by writing: $$ \bar{\alpha}_k = \cos\left(\frac{\theta}{2}\right)$$ $$ {\alpha}_k = e^{-i\phi}\sin\left(\frac{\theta}{2}\right)$$
We repeat this protocol iteratively for $k=0,1,...,2^N-1$.
