In Matlab, I simulated and graphed out the FFT of a signal composed of three pure sinusoids of distinct magnitudes and phases. Figuring out the scaling of the frequency axis is fairly straightforward - the number of bins in the frequency domain is equivalent to the sample size in the time domain. So bin #31 corresponds to $ 31 \times F_s/N = 31 \times 1000/1500 \approx 20 \operatorname{Hz}$.
What I don't understand is the scaling of the y-axis (magnitude or the phase). Given the definition of DFT and the math derivation from FS to DFT, there has to be a scaling factor of $\frac{1}{NT_s}$ between the digital amplitude and the analog amplitude:
$$ C_k \,=\, \frac{1}{T_0} \int_{-\frac{T_0}{2}}^\frac{T_0}{2} x(t) \, e^{-j \omega_0 kt} \operatorname{dt} \,=\, \frac{1}{T_0} \int_{T_0} x(t) \, e^{-j \frac{2\pi}{T_0} kt} \operatorname{dt} \,=\, \frac{1}{NT_s} \sum_{n=0}^{N-1} x(nT_s) \, e^{-j \frac{2\pi}{NT_s} k nT_s} \,=\, \frac{1}{NT_s} \underbrace{\sum_{n=0}^{N-1} x(nT_s) \, e^{-j \frac{2\pi}{N} k n}}_{\text{DFT}} $$
Take bin #31 as an example, the analog amplitude that corresponds to $ 20\operatorname{Hz} $ is $3$ (or two spectral lines each with height $1.5$ at $ \pm 20 \operatorname{Hz} $ because $ \cos(\theta) = 0.5e^{j\theta}+0.5^{-j\theta} $). If I scale the unscaled magnitude of $x[31]$ I get $ 2250 \times \frac{1}{NT_s} = 2250 \times \frac{1}{1500\times 1000} = 0.0015 $, which is obviously wrong. If I want to get to the value of $ 1.5 $, I have to divide the digital magnitude by $ 1500 $, which is exactly my sampling rate .... why?
Why do I rescale the magnitude by $ 1/N $ whereas the math derivation suggests a scaling factor of $ 1/(NT_s) $?
If I have to scale the magnitude by $ 1/N $ I would have thought the phase will need the same scaling - after all, both the magnitude and phase originate from the same complex spectrum equation. But to my surprise, the FFT gives me the exact phase value without any rescaling - at bin #31, the phase is $0.2$! How is that possible?
Matlab code:
fs=1000; % sampling @1000Hz
t=0:1/fs:1.5-1/fs; % from 0-1.499s in 0.001s interval so we have 1500 samples
f1=20;
f2=30;
f3=40;
x1=3*cos(2*pi*f1*t+0.2); % @f1=20Hz, amplitude=3, phase=0.2rad
x2=1*cos(2*pi*f2*t-0.3); % @f2=30Hz, amplitude=1, phase=-0.3rad
x3=2*cos(2*pi*f3*t+2.4); % @f3=40Hz, amplitude=2, phase=2.4rad
x=x1+x2+x3; % x has 1500 samples; x is a real number
X=fft(x); % X has 1500 "bins" of frequencies
X_mag = abs(X); % get magnitude spectrum from the complex spectrum; X_mag is a real number
figure(1) % unscaled magnitude spectrum
plot(X_mag);
figure(2) % scaled magnitude spectrum
plot(X_mag/1500);
figure(3) % unscaled phase spectrum
plot(angle(X));
FFT spectrum graphs: