SignalMixin.smooth_frames: use minimal data type for result.

Instead of always returning the result as an int64 or float64 array,
select the output type based on the types of the input arrays.

The output type should be the smallest type that has the correct
"kind" and is able to represent all input values.  For example, in
digital mode, if the input includes some int8 arrays and some int16
arrays, the result should be an int16 array.  In physical mode, if the
inputs are all float32, then the result will be float32; otherwise the
result will be float64.

However, although the output type should generally match the input
type, intermediate results may need to be stored as a different type.
For example, if the input and output are both int16, and one or more
signals have spf > 1 and use the entire 16-bit range, then the sum of
N samples will overflow an int16.  Previously, it was fine simply to
store the intermediate results in the output array itself, because the
output array was 64-bit, and no WFDB format has more than 32-bit
precision, and spf is (in practice) limited to at most 2**31-1.

For simplicity, continue using int64 or float64 as the intermediate
type, regardless of the actual input types and spf.

At the same time, we can also optimize the calculation slightly by
reshaping the input array and using np.sum, avoiding another Python
loop.

Author: Benjamin Moody

wfdb/io/_signal.py

@@ -830,32 +830,53 @@ def smooth_frames(self, sigtype='physical'):
         # Total samples per frame
         tspf = sum(spf)
 
+        # The output data type should be the smallest type that can
+        # represent any input sample value.  The intermediate data type
+        # must be able to represent the sum of spf[ch] sample values.
+
         if sigtype == 'physical':
             expanded_signal = self.e_p_signal
-            output_dtype = np.dtype('float64')
+            intermediate_dtype = np.dtype('float64')
+            allowed_dtypes = [
+                np.dtype('float32'),
+                np.dtype('float64'),
+            ]
         elif sigtype == 'digital':
             expanded_signal = self.e_d_signal
-            output_dtype = np.dtype('int64')
+            intermediate_dtype = np.dtype('int64')
+            allowed_dtypes = [
+                np.dtype('int8'),
+                np.dtype('int16'),
+                np.dtype('int32'),
+                np.dtype('int64'),
+            ]
         else:
             raise ValueError("sigtype must be 'physical' or 'digital'")
 
         n_sig = len(expanded_signal)
         sig_len = int(len(expanded_signal[0])/spf[0])
+        input_dtypes = set()
         for ch in range(n_sig):
             if len(expanded_signal[ch]) != sig_len * spf[ch]:
                 raise ValueError("length mismatch: signal %d has %d samples,"
                                  " expected %dx%d"
                                  % (ch, len(expanded_signal),
                                     sig_len, spf[ch]))
+            input_dtypes.add(expanded_signal[ch].dtype)
+
+        for output_dtype in allowed_dtypes:
+            if all(dt <= output_dtype for dt in input_dtypes):
+                break
+
         signal = np.zeros((sig_len, n_sig), dtype=output_dtype)
 
         for ch in range(n_sig):
             if spf[ch] == 1:
                 signal[:, ch] = expanded_signal[ch]
             else:
-                for frame in range(spf[ch]):
-                    signal[:, ch] += expanded_signal[ch][frame::spf[ch]]
-                signal[:, ch] = signal[:, ch] / spf[ch]
+                frames = expanded_signal[ch].reshape(-1, spf[ch])
+                signal_sum = np.sum(frames, axis=1, dtype=intermediate_dtype)
+                signal[:, ch] = signal_sum / spf[ch]
 
         return signal