How to add complex WGN to complex damped exponential and compute SNR?

I am new to the Signal processing domain, So my question might sound basic. I am trying to reproduce the results from the paper Decimative subspace-based parameter estimation techniques by Geert Morren, Philippe Lemmerling, Sabine Van Huffel.

I would like to first know what is the unit of the model given here, where $$a_k$$ amplitudes, $$\phi_k$$ phases, $$\alpha_k$$ damping factors, $$f_k$$ - frequencies in Hz and $$f_{sample}$$ - sampling frequency. $$\begin{equation} x(t) = \sum_{k=1}^{K} (a_k e^{(j\phi_k)})(e^{\{(j2\pi f_k - \alpha_k)/f_{sample}\}t}) + Noise(t), \quad t = 0,1,...,N-1 \end{equation}$$

I am adding circularly symmetric WGN to $$x(t)$$ using the wgn function from MATLAB (shown below in code), which is in dBW unit. I am not sure whether this is correct, and I am computing SNR also using a snr MATLAB function. However, I want to use the variance of noise and scaling of noise to determine SNR. Could someone explain me clearly how to do this?

f       = [0.02, 0.0205];
a       = [1, 1];
alpha   = [0, 0];
phi     = [0, 0];
K       = length(f);
del_t   = 1;
N       = 1000;
t       = (0:del_t:N-1)';
size_t  = length(t);
f_sample= 1/del_t;
Nruns   = 40;
x       = zeros(size_t,1);
SNR     = zeros(1,Nruns);

for iter = 1:Nruns
noise = wgn(size_t,1,randi(20),'complex');
for n = drange(1:size_t)
for k = 1:K
x(n) = x(n) + a(k)*exp(phi(k)*1j)*(exp((2*pi*f(k)*1j - alpha(k))/f_sample)).^(n-1);
end
x(n) = x(n) + noise(n);
end
SNR(iter) = abs(snr((x-noise),noise));
end

One suggestion is to leverage the vector processing in Matlab and to not add the noise in a for loop, but given the vector of your signal of length $$N$$ you can simply create the white Gaussian Noise using:

noise = std * sqrt(2)/2 * randn(N,1) + std * sqrt(2)/2* 1j*randn(N,1)

(or randn(1,N) if you prefer row vectors instead of column vectors)

Where std represents the total standard deviation of the noise.

Next be aware that the total power as a white noise and the square of the standard deviation std in whatever units are chosen is distributed evenly over the entire frequency range given by -$$f_s/2$$ to +$$f_s/2$$ where $$f_s$$ is the sampling rate. SNR of interest is typically that portion of the noise that shares the same spectrum as the signal. (So if the signal only occupied 1/4 of the total spectrum as sampled, then the noise would be 1/4 or 6 dB lower than that predicted by the total variance of noise).