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.
Then the noise can be added to your signal simply using
y = x + 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).