Calculation of motion RAOs (Response Amplitude Operators) using Fourier Transform (MATLAB)

An example of using MATLAB to derive the motion RAOs of a body induced by a regular wave. The method can also be applied for obtaining the dominant component in a signal.

 
load ex_fromPointTracker20190416.mat % data from image processing
Aamp=[1.1859 1.3743 1.3733 1.6982 1.4289 1.4463 1.6919]; % wave cond of no WEC no wind
Prd=[0.70 0.900 1.00 1.2 1.40 1.5 1.8];
omega=2*pi./Prd;
g=9.81;
k=omega.^2/g; % wave number
den=k.*Aamp*1e-2; % for dimensionless
L1=-0.2; 
ratio=L1/(exCase{1,1}.pos(1,2)-exCase{1,2}.pos(1,2))*100;
L2=70; % unit: cm
f=59.94;Tc=1/f;
for i=1:5
    pitch=atan2(-(exCase{i,2}.pos(:,2)-exCase{i,1}.pos(:,2)),exCase{i,2}.pos(:,1)-exCase{i,1}.pos(:,1));
    pitch=pitch-pi/2;
    pitch=detrend(pitch);
    t=0:Tc:(length(pitch)-1)*Tc;
    figure(3);plot(t,pitch);
    heave=(exCase{i,1}.pos(:,2))*ratio./Aamp(i);
    figure(2);plot(t,heave);
    surge=((exCase{i,1}.pos(:,1)-exCase{i,1}.pos(1,1))*ratio-L2*sin(-pitch))./Aamp(i);
    surge=detrend(surge);
    heave=detrend(heave);
    figure(1);plot(t,surge);
    [f,amp]=fft_freq(surge,Tc);
    figure(4);plot(f,amp);
    surgeRao(i)=interp1(f,amp,1/Prd(i));
    [f,amp]=fft_freq(heave,Tc);
    figure(5);plot(f,amp);
    heaveRao(i)=interp1(f,amp,1/Prd(i));
    pitch=pitch*180./pi/den(i);
    [f,amp]=fft_freq(pitch,Tc);
    figure(6);plot(f,amp);
    pitchRao(i)=interp1(f,amp,1/Prd(i));
end
    

casex=1:5;
figure(11);
plot(casex,surgeRao);
figure(12);
plot(casex,heaveRao);
figure(13);
plot(casex,pitchRao);

Codes can be downloaded from motionRAO.zip