function [yt] = spectrum2tseries(f,P,t)


%Units of P: Power/Hz, where power is amplitude squared
%Units of f: Hz


df = f(2)-f(1);


zs = (2*P*df).^0.5; %spectral energy to eta (wave amplitude)

g = 9.81; 
k = (2*pi*f).^2/g; %wavenumber from Dispersion relation

%to generate random numbers in the range (0,2*pi) 
n_phases = length(f);
phasea = 0;
phaseb = 2*pi;
         
phases = phasea + (phaseb-phasea).*rand(n_phases,1);

  x = 0;

for j = 1:length(t);
     yt(j) = zs(1)*cos(k(1)*x - 2*pi*f(1)*t(j) + phases(1));
  
     for i = 2:length(f)
          yt(j) = yt(j) +  zs(i)*cos(k(i)*x - 2*pi*f(i)*t(j) + phases(i));
     end
     
end