Downsample a 1D numpy array

Question

I have a 1-d numpy array which I would like to downsample. Any of the following methods are acceptable if the downsampling raster doesn't perfectly fit the data:

• overlap downsample intervals
• convert whatever number of values remains at the end to a separate downsampled value
• interpolate to fit raster

basically if I have

1 2 6 2 1

and I am downsampling by a factor of 3, all of the following are ok:

3 3

3 1.5

or whatever an interpolation would give me here.

I'm just looking for the fastest/easiest way to do this.

I found scipy.signal.decimate, but that sounds like it decimates the values (takes them out as needed and only leaves one in X). scipy.signal.resample seems to have the right name, but I do not understand where they are going with the whole fourier thing in the description. My signal is not particularly periodic.

Could you give me a hand here? This seems like a really simple task to do, but all these functions are quite intricate...

Solution 1

In the simple case where your array's size is divisible by the downsampling factor (R), you can reshape your array, and take the mean along the new axis:

import numpy as np
a = np.array([1.,2,6,2,1,7])
R = 3
a.reshape(-1, R)
=> array([[ 1.,  2.,  6.],
         [ 2.,  1.,  7.]])

a.reshape(-1, R).mean(axis=1)
=> array([ 3.        ,  3.33333333])

In the general case, you can pad your array with NaNs to a size divisible by R, and take the mean using scipy.nanmean.

import math, scipy
b = np.append(a, [ 4 ])
b.shape
=> (7,)
pad_size = math.ceil(float(b.size)/R)*R - b.size
b_padded = np.append(b, np.zeros(pad_size)*np.NaN)
b_padded.shape
=> (9,)
scipy.nanmean(b_padded.reshape(-1,R), axis=1)
=> array([ 3.        ,  3.33333333,  4.])

Solution 2

Here are a few approaches using either linear interpolation or the Fourier method. These methods support upsampling as well as downsampling.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample
from scipy.interpolate import interp1d

def ResampleLinear1D(original, targetLen):
    original = np.array(original, dtype=np.float)
    index_arr = np.linspace(0, len(original)-1, num=targetLen, dtype=np.float)
    index_floor = np.array(index_arr, dtype=np.int) #Round down
    index_ceil = index_floor + 1
    index_rem = index_arr - index_floor #Remain

    val1 = original[index_floor]
    val2 = original[index_ceil % len(original)]
    interp = val1 * (1.0-index_rem) + val2 * index_rem
    assert(len(interp) == targetLen)
    return interp

if __name__=="__main__":

    original = np.sin(np.arange(256)/10.0)
    targetLen = 100

    # Method 1: Use scipy interp1d (linear interpolation)
    # This is the simplest conceptually as it just uses linear interpolation. Scipy
    # also offers a range of other interpolation methods.
    f = interp1d(np.arange(256), original, 'linear')
    plt.plot(np.apply_along_axis(f, 0, np.linspace(0, 255, num=targetLen)))

    # Method 2: Use numpy to do linear interpolation
    # If you don't have scipy, you can do it in numpy with the above function
    plt.plot(ResampleLinear1D(original, targetLen))

    # Method 3: Use scipy's resample
    # Converts the signal to frequency space (Fourier method), then back. This
    # works efficiently on periodic functions but poorly on non-periodic functions.
    plt.plot(resample(original, targetLen))

    plt.show()