Generating irregular waves in MATLAB

So, I've been trying to reproduce in MATLAB the strategy for generating irregular waves from this tech report: Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine

But I am missing something: I have to define a bin size in which to subdivide the angular frequencies used to gather values from the PSD and generate the white gaussian noise, but it should not affect the end result. The issue is, it was affecting it by changing the overall amplitude and frequency of the waves.

So I went back to studying and came up with solutions to both issues:

• Dividing the PSD by the bin size (step) for the amplitude
• Having the period of my simulation changing proportionally with the bin size (step).

I mean, I can see why sampling more points from the PSD leads to more data when in time domain. But I could not find a reasonable explanation to why dividing the PSD by the bin size yields the same amplitude when you change it. Bellow follows my code, with 'noise' and 'spectra' being the functions I used to create the arrays for the white gaussian noise and the PSD, respectivelly.

%Determining the frequency range
lim = 1.2;
step = 0.01;
w = step:step:lim;

%Preparing random variates for the WGN generation
rng('shuffle');
u1 = rand(1,size(w,2));
rng('shuffle');
u2 = rand(1,size(w,2));

%Calling the noise fun. for every w and mirroring it
n = noise(u1,u2);
rn = fliplr(conj(n));
wgn = [rn, 0, n];

%Calling the spectra fun. for every w and mirroring it
s = spectra(w);
rs = rot90(s,2);
sout = [rs,0,s];

%2-sided freq. and time definition
w = -lim:step:lim;
Ts = 2*pi/lim;
tlim = Ts*length(sout)-Ts;
t = 0:Ts:tlim;

%Processing the data, taking it to time domain
pre_waves = wgn.*sqrt(pi.*sout./step);
waves = ifft(ifftshift(pre_waves), 'symmetric');

%Plotting
figure(3)
plot(t,waves)
grid