% Simple simulation of a complex Gaussian random process

N = 512;                % Number of temporal samples to generate
NLags = 512;            % Numer of lags in the corrlation function
SNR = 20;               % Signal-to-noise ratio
v = 0.20;               % Radial velocity of scatterer
PRF = 4e3;              % Pulse repitition frequency
f0 = 5e6;               % Pulse center frequency
Np = 4.5;               % Number of pulse periods
c = 1540;               % Speed of sound

fd = (2*f0*v/c)/PRF;    % Theoretical Doppler shift
L = Np*c/f0/2;          % Pulse length
Tr = L/v;               % Transit time radially
Tprf = 1/PRF;
s2 = Tr/Tprf;           % Normalized transit time

m = 0:NLags-1;
R = 10^(SNR/10)*exp(-0.5*(m/s2).^2).*exp(i*2*pi*fd*m); % Gaussian correlation function given by the ultrasound parameters

% Plot correlation function
figure(102);
plot(m,abs(R)/max(abs(R)));
xlim([0 100]);
title('Correlation function');

w = (freqspace(N,'full')-1)*pi;

% Fourier transform correlation function
G = fft(R,N).';
% Generate white noise
Z = randn(N,1) + i*randn(N,1);
% Shape white noise
Z = sqrt(real(G)).*Z;
% Invers transform to time domain
z = ifft(Z);
% Estimate power spectrum from simulated signal
Gest = fftshift(abs(fft(z,N)).^2);

% Compare true versus simulated spectrum
figure(100);
clf;
hold on;
plot(w,10*log10(fftshift(real(G))),'LineWidth',1.5);
plot(w,10*log10(Gest),'r');
hold off;
xlim([-pi pi]);
%ylim([-80 40]);