Compute surface elevation using IFFT

.gitignore

@@ -32,3 +32,6 @@ $ cat .gitignore
 # Exemptions
 !**/examples/data/wave/*.mat
 !**/tests/data/wave/*.mat
+
+# Files created during tests
+mhkit/tests/wave/plots/

mhkit/tests/wave/test_resource_spectrum.py

@@ -32,12 +32,11 @@ class TestResourceSpectrum(unittest.TestCase):
 
     @classmethod
     def setUpClass(self):
-        omega = np.arange(0.1,3.5,0.01)
-        self.f = omega/(2*np.pi)
+        Trep = 600
+        df = 1 / Trep
+        self.f = np.arange(0, 1, df)
         self.Hs = 2.5
         self.Tp = 8
-        df = self.f[1] - self.f[0]
-        Trep = 1/df
         self.t = np.arange(0, Trep, 0.05)
 
     @classmethod
@@ -55,6 +54,17 @@ def test_pierson_moskowitz_spectrum(self):
         self.assertLess(errorHm0, 0.01)
         self.assertLess(errorTp0, 0.01)
 
+    def test_pierson_moskowitz_spectrum_zero_freq(self):
+        df = 0.1
+        f_zero = np.arange(0, 1, df)
+        f_nonzero = np.arange(df, 1, df)
+
+        S_zero = wave.resource.pierson_moskowitz_spectrum(f_zero, self.Tp, self.Hs)
+        S_nonzero = wave.resource.pierson_moskowitz_spectrum(f_nonzero, self.Tp, self.Hs)
+
+        self.assertEqual(S_zero.values.squeeze()[0], 0.0)
+        self.assertGreater(S_nonzero.values.squeeze()[0], 0.0)
+
     def test_jonswap_spectrum(self):
         S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
         Hm0 = wave.resource.significant_wave_height(S).iloc[0,0]
@@ -66,6 +76,17 @@ def test_jonswap_spectrum(self):
         self.assertLess(errorHm0, 0.01)
         self.assertLess(errorTp0, 0.01)
 
+    def test_jonswap_spectrum_zero_freq(self):
+        df = 0.1
+        f_zero = np.arange(0, 1, df)
+        f_nonzero = np.arange(df, 1, df)
+
+        S_zero = wave.resource.jonswap_spectrum(f_zero, self.Tp, self.Hs)
+        S_nonzero = wave.resource.jonswap_spectrum(f_nonzero, self.Tp, self.Hs)
+
+        self.assertEqual(S_zero.values.squeeze()[0], 0.0)
+        self.assertGreater(S_nonzero.values.squeeze()[0], 0.0)
+
     def test_surface_elevation_phases_np_and_pd(self):
         S0 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
         S1 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs*1.1)
@@ -129,6 +150,14 @@ def test_surface_elevation_rmse(self):
 
         self.assertLess(rmse_sum, 0.02)
 
+    def test_ifft_sum_of_sines(self):
+        S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
+
+        eta_ifft = wave.resource.surface_elevation(S, self.t, seed=1, method='ifft')
+        eta_sos = wave.resource.surface_elevation(S, self.t, seed=1, method='sum_of_sines')
+
+        assert_allclose(eta_ifft, eta_sos)     
+
     def test_plot_spectrum(self):
         filename = abspath(join(plotdir, 'wave_plot_spectrum.png'))
         if isfile(filename):

mhkit/wave/resource.py

@@ -90,7 +90,17 @@ def pierson_moskowitz_spectrum(f, Tp, Hs):
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
     A_PM = B_PM*(Hs/2)**2
-    Sf  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
+
+    # Avoid a divide by zero if the 0 frequency is provided
+    # The zero frequency should always have 0 amplitude, otherwise
+    # we end up with a mean offset when computing the surface elevation.
+    Sf = np.zeros(f.size)
+    if f[0] == 0.0:
+        inds = range(1, f.size)
+    else:
+        inds = range(0, f.size)
+
+    Sf[inds]  = A_PM*f[inds]**(-5)*np.exp(-B_PM*f[inds]**(-4))
 
     col_name = 'Pierson-Moskowitz ('+str(Tp)+'s)'
     S = pd.DataFrame(Sf, index=f, columns=[col_name])
@@ -132,7 +142,17 @@ def jonswap_spectrum(f, Tp, Hs, gamma=None):
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
     A_PM = B_PM*(Hs/2)**2
-    S_f  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
+
+    # Avoid a divide by zero if the 0 frequency is provided
+    # The zero frequency should always have 0 amplitude, otherwise
+    # we end up with a mean offset when computing the surface elevation.
+    S_f = np.zeros(f.size)
+    if f[0] == 0.0:
+        inds = range(1, f.size)
+    else:
+        inds = range(0, f.size)
+
+    S_f[inds]  = A_PM*f[inds]**(-5)*np.exp(-B_PM*f[inds]**(-4))
 
     if not gamma:
         TpsqrtHs = Tp/np.sqrt(Hs);
@@ -162,7 +182,7 @@ def jonswap_spectrum(f, Tp, Hs, gamma=None):
     return S
 
 ### Metrics
-def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None):
+def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None, method='ifft'):
     """
     Calculates wave elevation time-series from spectrum
 
@@ -175,9 +195,18 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
         for example, time = np.arange(0,100,0.01)
     seed: int (optional)
         Random seed
+    frequency_bins: numpy array or pandas DataFrame (optional)
+        Bin widths for frequency of S. Required for unevenly sized bins
     phases: numpy array or pandas DataFrame (optional)
         Explicit phases for frequency components (overrides seed)
         for example, phases = np.random.rand(len(S)) * 2 * np.pi
+    method: str (optional)
+        Method used to calculate the surface elevation. 'ifft'
+        (Inverse Fast Fourier Transform) used by default if the
+        given frequency_bins==None.
+        'sum_of_sines' explicitly sums each frequency component
+        and used by default if frequency_bins are provided.
+        The 'ifft' method is significantly faster.
 
     Returns
     ---------
@@ -194,13 +223,21 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
             "frequency_bins must be of type None, np.ndarray, or pd,DataFrame")
     assert isinstance(phases, (type(None), np.ndarray, pd.DataFrame)), (
             'phases must be of type None, np.ndarray, or pd,DataFrame')
+    assert isinstance(method, str)
 
     if frequency_bins is not None:
         assert frequency_bins.squeeze().shape == (S.squeeze().shape[0],),(
             'shape of frequency_bins must match shape of S')
     if phases is not None:
         assert phases.squeeze().shape == S.squeeze().shape,(
             'shape of phases must match shape of S')
+
+    if method is not None:
+        assert method == 'ifft' or method == 'sum_of_sines',(
+            f"unknown method {method}, options are 'ifft' or 'sum_of_sines'")
+
+    if method == 'ifft':
+        assert S.index.values[0] == 0, ('ifft method must have zero frequency defined')        
 
     f = pd.Series(S.index)
     f.index = f
@@ -209,10 +246,12 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
         assert np.allclose(f.diff()[1:], delta_f)
     elif isinstance(frequency_bins, np.ndarray):
         delta_f = pd.Series(frequency_bins, index=S.index)
+        method = 'sum_of_sines'
     elif isinstance(frequency_bins, pd.DataFrame):
         assert len(frequency_bins.columns) == 1, ('frequency_bins must only'
                 'contain 1 column')
         delta_f = frequency_bins.squeeze()
+        method = 'sum_of_sines'
 
     if phases is None:
         np.random.seed(seed)
@@ -231,17 +270,28 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
     A = A.multiply(delta_f, axis=0)
     A = np.sqrt(A)
 
-    # Product of omega and time
-    B = np.outer(time_index, omega)
-    B = B.reshape((len(time_index), len(omega)))
-    B = pd.DataFrame(B, index=time_index, columns=omega.index)
-
-    # wave elevation
-    eta = pd.DataFrame(columns=S.columns, index=time_index)
-    for mcol in eta.columns:
-        C = np.cos(B+phase[mcol])
-        C = pd.DataFrame(C, index=time_index, columns=omega.index)
-        eta[mcol] = (C*A[mcol]).sum(axis=1)
+    if method == 'ifft':
+        A_cmplx = A * (np.cos(phase) + 1j*np.sin(phase))
+
+        def func(v):
+            eta = np.fft.irfft(0.5 * v.values.squeeze() * time_index.size, time_index.size)
+            return pd.Series(data=eta, index=time_index)
+
+        eta = A_cmplx.apply(func)
+
+    elif method == 'sum_of_sines':
+        # Product of omega and time
+        B = np.outer(time_index, omega)
+        B = B.reshape((len(time_index), len(omega)))
+        B = pd.DataFrame(B, index=time_index, columns=omega.index)
+
+        # wave elevation
+        eta = pd.DataFrame(columns=S.columns, index=time_index)
+        for mcol in eta.columns:
+            C = np.cos(B+phase[mcol])
+            C = pd.DataFrame(C, index=time_index, columns=omega.index)
+            eta[mcol] = (C*A[mcol]).sum(axis=1)
+
     return eta