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: Select all

`%Read the WAMIT output file for the Barge Platform Wave-Excitation Force`

BX = xlsread('C:\Users\Path\To\Barge.3\File\In\xlsx\Format','A1:G22201');

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 = 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');

pre_force5 = wgn.*sqrt(pi.*sout./step).*X5;

force5 = ifft(pre_force5,'symmetric');

%Plotting

figure(3)

plot(t,waves,'b', t, force5, 'r')

grid