rel_calib_stack: fix bug when using non-default windowing overlap fraction

CHANGELOG.txt

@@ -142,6 +142,8 @@ master: (doi: 10.5281/zenodo.165135)
      alongside gaps (see #1366)
    * obspy-plot now has option to disable min/max plot (see #1583)
  - obspy.signal:
+   * fixed a bug in calibration.rel_calib_stack (resulting amplitude response
+     had wrong scaling if using non-default "overlap_fraction", see #1821)
    * New obspy.signal.quality_control module to compute quality metrics from
      MiniSEED files. (see #1141)
    * New correlate function for calculating the cross-correlation function

obspy/signal/calibration.py

@@ -73,7 +73,6 @@ def rel_calib_stack(st1, st2, calib_file, window_len, overlap_frac=0.5,
         msg = "Traces don't have the same sampling rate!"
         raise ValueError(msg)
     else:
-        ndat1 = st1[0].stats.npts
         sampfreq = st1[0].stats.sampling_rate
 
     # read waveforms
@@ -88,14 +87,16 @@ def rel_calib_stack(st1, st2, calib_file, window_len, overlap_frac=0.5,
     gg, _freq = _calc_resp(calib_file, nfft, sampfreq)
 
     # calculate number of windows and overlap
-    nwin = int(np.floor((ndat1 - nfft) / (nfft / 2)) + 1)
     noverlap = nfft * overlap_frac
 
     auto, _freq, _t = \
         spectral_helper(tr1, tr1, NFFT=nfft, Fs=sampfreq, noverlap=noverlap)
     cross, freq, _t = \
         spectral_helper(tr2, tr1, NFFT=nfft, Fs=sampfreq, noverlap=noverlap)
 
+    # get number of windows that were actually computed inside FFT routine
+    nwin = auto.shape[1]
+
     res = (cross / auto).sum(axis=1) * gg
 
     # The first item might be zero. Problems with phase calculations.

obspy/signal/tests/test_calibration.py

@@ -89,6 +89,35 @@ def test_relcal_using_traces(self):
         np.testing.assert_array_almost_equal(amp, amp2, decimal=4)
         np.testing.assert_array_almost_equal(phase, phase2, decimal=4)
 
+    def test_relcal_different_overlaps(self):
+        """
+        Tests using different window overlap percentages.
+
+        Regression test for bug #1821.
+        """
+        st1 = read(os.path.join(self.path, 'ref_STS2'))
+        st2 = read(os.path.join(self.path, 'ref_unknown'))
+        calfile = os.path.join(self.path, 'STS2_simp.cal')
+
+        def median_amplitude_plateau(freq, amp):
+            # resulting response is pretty much flat in this frequency range
+            return np.median(amp[(freq >= 0.3) & (freq <= 3)])
+
+        # correct results using default overlap fraction of 0.5
+        freq, amp, phase = rel_calib_stack(
+            st1, st2, calfile, 20, smooth=10, overlap_frac=0.5,
+            save_data=False)
+        amp_expected = median_amplitude_plateau(freq, amp)
+        for overlap in np.linspace(0.1, 0.9, 5):
+            freq2, amp2, phase2 = rel_calib_stack(
+                st1, st2, calfile, 20, smooth=10, overlap_frac=overlap,
+                save_data=False)
+            amp_got = median_amplitude_plateau(freq2, amp2)
+            percentual_difference = abs(
+                (amp_expected - amp_got) / amp_expected)
+            # make sure results are close for any overlap choice
+            self.assertTrue(percentual_difference < 0.01)
+
     def test_relcal_using_dict(self):
         """
         Tests using paz dictionary instead of a gse2 file.