Magnitude of Power Spectral Density with Different Sampling Frequency

Question

For exponential signals (sine or cosine), if the sampling frequency $f_s$ is equal to the length of signal, $N$, the magnitude of psd for each sine signal is proportional to amplitude in time domain, while if $f_s$ is not equal to length of signal $N$, for example here, $f_s=1000 \text{Hz}$, the magnitude of sine signals in psd (frequency domain) is not proportional to their original signal. The MATLAB code is provided.

Can you tell me the reason? Is this something related with leakage?

N=256;
x0=zeros(N,1);
fs=1000;
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X).^2)

Answer 1

The issue (beside the scaling) here is that your selected frequencies do not align with your FFT grid and that the energy of each complex exponential is smeared out over multiple FFT bins. Your frequency resolution is given by sample rate divide by FFT length. In your case that's 3.906 Hz. Only for frequencies that are an integer multiple of this will all energy fall into a single FFT bin. An easy way to achieve this is to bump the sample rate up to 1024 Hz. Then your frequency resolution is exactly 4 Hz and all your frequencies are an integer multiple of that. Once you fix the scaling as well, it works fine like this:

N=256;
x0=zeros(N,1);
fs=1024; % Fix 1: get to exactly 4Hz frequency resolution
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X)/N); % Fix 2: absolute value only, rescaled

Answer 2

I am not sure I understand your question, nevertheless, I would like to point out that matlab computes fft considering this expression: $$ D_{n}=\sum_{k=0}^{K-1}x(k)\cdot e^{-j\cdot n\cdot2\cdot\pi\cdot\frac{k}{K}} $$

Which as you know it misses a factor of $$\frac{1}{K} $$

That is: $$ D_{n}=\frac{1}{K}\sum_{k=0}^{K-1}x(k)\cdot e^{-j\cdot n\cdot2\cdot\pi\cdot\frac{k}{K}} $$

So, in your matlab code you should call: X=fft(x0)/N;

This book is great addressing this topic. Go to section 2.10 (Numerical Computation of D_n) on page 60, you will also find there some useful matlab code.