WVM issue in scaling.py

pvlib/scaling.py

@@ -159,7 +159,12 @@ def latlon_to_xy(coordinates):
 
 def _compute_wavelet(clearsky_index, dt=None):
     """
-    Compute the wavelet transform on the input clear_sky time series.
+    Compute the wavelet transform on the input clear_sky time series. Uses a
+    Haar wavelet based on moving average. Output resembles a Stationary Wavelet
+    Transform (SWT), but has shifted time positions on the averages relative to
+    output from other wavelet codes. Returns one level of approximation
+    coefficient (CAn) and n levels of detail coefficients (CD1, CD2, ...,
+    CDn-1, CDn). Maximum scale is 4096 s (12 levels at a signal dt of 1s).
 
     Parameters
     ----------
@@ -174,7 +179,8 @@ def _compute_wavelet(clearsky_index, dt=None):
     Returns
     -------
     wavelet: numeric
-        The individual wavelets for the time series
+        The individual wavelets for the time series. Format follows increasing
+        scale (decreasing frequency): [CD1, CD2, ..., CDn, CAn]
 
     tmscales: numeric
         The timescales associated with the wavelets in seconds [s]
@@ -211,29 +217,31 @@ def _compute_wavelet(clearsky_index, dt=None):
     min_tmscale = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet timescale
     max_tmscale = int(12 - min_tmscale)  # maximum wavelet timescale
 
-    tmscales = np.zeros(max_tmscale)
-    csi_mean = np.zeros([max_tmscale, len(cs_long)])
-    # Loop for all time scales we will consider
-    for i in np.arange(0, max_tmscale):
-        j = i+1
-        tmscales[i] = 2**j * dt  # Wavelet integration time scale
-        intvlen = 2**j  # Wavelet integration time series interval
+    # Compute approximation coefficients via moving average (Haar wavelet)
+    csi_mean = np.zeros([max_tmscale+1, len(cs_long)])
+    csi_mean[0, :] = cs_long.values.flatten()  # CA0, the original signal
+    for i in np.arange(1, max_tmscale+1):
+        intvlen = 2**i  # Wavelet integration time series interval
         # Rolling average, retains only lower frequencies than interval
         df = cs_long.rolling(window=intvlen, center=True, min_periods=1).mean()
         # Fill nan's in both directions
         df = df.fillna(method='bfill').fillna(method='ffill')
         # Pop values back out of the dataframe and store
-        csi_mean[i, :] = df.values.flatten()
+        csi_mean[i, :] = df.values.flatten()  # Approximation Coeffs CA0,CA1,...
 
-    # Calculate the wavelets by isolating the rolling mean frequency ranges
+    # Calculate detail coeffs by difference between approximation levels
     wavelet_long = np.zeros(csi_mean.shape)
-    for i in np.arange(0, max_tmscale-1):
+    tmscales = np.zeros(max_tmscale+1)
+    for i in np.arange(0, max_tmscale):
+        # Compute detail coeffs, beginning with CD1 through CDn
         wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :]
-    wavelet_long[max_tmscale-1, :] = csi_mean[max_tmscale-1, :]  # Lowest freq
+        tmscales[i] = 2**(i+1) * dt  # Wavelet integration time scale
+    wavelet_long[max_tmscale, :] = csi_mean[max_tmscale, :]  # Lowest freq (CAn)
+    tmscales[max_tmscale] = tmscales[max_tmscale-1]  # CAn and CDn share scale
 
     # Clip off the padding and just return the original time window
-    wavelet = np.zeros([max_tmscale, len(vals)])
-    for i in np.arange(0, max_tmscale):
-        wavelet[i, :] = wavelet_long[i, len(vals)+1: 2*len(vals)+1]
+    wavelet = np.zeros([max_tmscale+1, len(vals)])
+    for i in np.arange(0, max_tmscale+1):
+        wavelet[i, :] = wavelet_long[i, len(vals): 2*len(vals)]
 
     return wavelet, tmscales

pvlib/tests/test_scaling.py

@@ -48,21 +48,21 @@ def positions():
 @pytest.fixture
 def expect_tmscale():
     # Expected timescales for dt = 1
-    return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+    return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 4096]
 
 
 @pytest.fixture
 def expect_wavelet():
-    # Expected wavelet for indices 5000:5004 for clear_sky_index above (Matlab)
-    return np.array([[-0.025, 0.05, 0., -0.05, 0.025],
-                     [0.025, 0., 0., 0., -0.025],
-                     [0., 0., 0., 0., 0.]])
+    # Expected wavelet for indices 5000:5004 for clear_sky_index above (Manual)
+    return np.array([[0., 0., 0.05, -0.1, 0.05],
+                     [0., -0.025, 0.05, 0., -0.05],
+                     [0., 0.025, 0., 0., 0.]])
 
 
 @pytest.fixture
 def expect_cs_smooth():
     # Expected smoothed clear sky index for indices 5000:5004 (Matlab)
-    return np.array([1., 1.0289, 1., 0.9711, 1.])
+    return np.array([1., 1., 1.0577, 0.9423, 1.])
 
 
 def test_latlon_to_xy_zero():
@@ -97,6 +97,14 @@ def test_compute_wavelet_series(clear_sky_index, time,
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
+def test_compute_wavelet_dwttheory(clear_sky_index, time,
+                                   expect_tmscale, expect_wavelet):
+    # Confirm detail coeffs sum to original signal
+    csi_series = pd.Series(clear_sky_index, index=time)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(np.sum(wavelet, 0), csi_series)
+
+
 def test_compute_wavelet_series_numindex(clear_sky_index, time,
                                          expect_tmscale, expect_wavelet):
     dtindex = pd.to_datetime(time, unit='s')