clc;clear;
close all;
Fs= 20;  % Fs is the sampling rate in one second (in Hz)

FilePath={'Test25.out'};
FilesDesc={};
ReferenceFile=[];
ShowLegend=[1];
CustomHdr={[],8,7,8};
PlotPSDs=[true];
OnePlot=[0];


[outData]=PlotFASToutput(FilePath,FilesDesc,ReferenceFile,{'Wave1Elev'},ShowLegend,CustomHdr,PlotPSDs,OnePlot)
n=length(outData{1,2}(:,1));
if mod(n,2)==1
   n=n-1;
end
[f1, pxx] = Jason_PSD( outData{1,2}(1:n,1), 1/(outData{1,1}(n) - outData{1,1}(1)) );  % calculate auto power density
x=outData{1,2}(1:n,1);


[outData]=PlotFASToutput(FilePath,FilesDesc,ReferenceFile,{'PtfmSurge'},ShowLegend,CustomHdr,PlotPSDs,OnePlot);
y=outData{1,2}(1:n,1);
pxy=cpsd(x,y,[],[],[],Fs); % cpsd is a subroutine self-contained in Matlab to calculate the cross power density

RAO=pxy/pxx;