# Irregular Waves and Pitch Momentum in MATLAB

Hello,

I’ ve been trying to replicate the strategy used to generate irregular waves and the diffraction pitch momentum as described in: https://www.nrel.gov/docs/fy08osti/41958.pdf using MATLAB, for academic purposes. However, I’m not sure what I am doing wrong, as it seems that, when I increase the sampling size for the angular frequency (the ‘step’ variable in my script), it affects both amplitude and frequency of the waves generated, when it (probably, I guess) shouldn’t happen.
What seems to work is dividing the PSD by the step size and making it so that the simulation time increases as the step size decreases. The second one makes sense to me, since in the inverse fast fourier transform, more points in the frequency domain equals more time before ‘repeating’ itself. But I could not find an explanation to the amplitude issue. Following is my code, and examples of graphs with and without corrections. Any help is apprecriated.

[code]%Read the WAMIT output file for the Barge Platform Wave-Excitation Force
B5 = BX(BX(:,2)==0 & BX(:,3)==5,:);
B5(:,1) = 2*pi./B5(:,1);

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

%Reading only data regarding Beta = 0 and Pitch motion, summing both real
%and imaginary parts
ReX5 = interp1(B5(:,1), B5(:,6), w, ‘linear’);
ImX5 = interp1(B5(:,1), B5(:,7), w, ‘linear’);
X5 = ReX5+1i*ImX5;

%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];

%Taking the wave-excitation vector and mirroring it
X5n = fliplr(conj(n));
X5 = [X5n, 0, X5];

%New boundaries defined
w = -lim:step:lim;
Ts = 2pi/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’);

pre_force5 = wgn.*sqrt(pi.*sout./step).*X5;
force5 = ifft(pre_force5,‘symmetric’);

%Plotting
figure(3)
plot(t,waves,‘b’, t, force5, ‘r’)
grid[/code]

Dear Leonardo,

I haven’t reviewed your MATLAB code in detail, and I don’t see any results attached to review. Regardless, in the FFT method, the frequency step is related to the simulation length and the time step is related to maximum frequency, i.e.:

domega = 2*pi/t_max
omega_max = pi/dt

(The frequencies here are in rad/s, not Hz.)

This is more clearly explained in my 2009 Wind Energy journal paper, which is based on the report you mention: onlinelibrary.wiley.com/doi/abs/10.1002/we.347.

When changing the simulation length, the frequency step should change, which means the random phases tied to each frequent component will change and you’ll get a different realization of the wave time series.

Best regards,

Dear Jason,

Thank you for your reply. I have implemented the relation between time and frequency you cited, and it seems to be working:

But the amplitude seems to be the issue, since I’ve used the same JONSWAP spectra from your paper, the significant height should look like:

In this simulation, I’ve divided the spectrum by the ‘domega’ frequency step, and it seems to work, but I don’t know why, as I’ve done this out of intuition. In the meantime, I’ll be reading the paper you suggested.

Dear Leanardo,

The DFT of the White Guassian Noise using the Box-Muller method is not stated in the Wind Energy paper I referenced above. But this is:

W[k] = SQRT(N/2)SQRT(-2LN(U1[k]))EXP(j2piU2[k])

where U1[k] and U2[k] are independent and uniformly distributed random variates (random numbers between zero and one). U2[k] randomizes the phase and U1[k] randomizes the amplitude. If the amplitudes are not to be randomized (WaveNDAmp = FALSE in HydroDyn), then SQRT(-2*LN(U1[k])) is replaced by its mean value, SQRT(2).

Best regards,