Oscillation when using scipy.signal.decimate with ftype='fir'

I am porting some code from Matlab (v2019b) to python (3.9.12, scipy 1.7.3). I am seeing differences between the results from the Matlab call decimate(A, subsample, 'fir') and scipy.signal.decimate(A, subsample, ftype='fir'). I cannot share the code I am working on, but here is a toy example:

Matlab:

artest = 1:109;
artestd = decimate(artest, 2, 'fir');
artestd(40:end)

ans =

Columns 1 through 12

79.0000 81.0000 83.0000 85.0000 87.0000 89.0000 91.0000 93.0000 95.0000 97.0000 99.0000 101.0000

Columns 13 through 16

103.0000 105.0000 107.0000 109.0000

Python:

import numpy as np
from scipy import signal as sg
artest = np.arange(1,110)
artestd = sg.decimate(artest, 2, ftype='fir')
artestd[39:]

Out[6]:

array([ 79.        ,  81.        ,  83.        ,  85.        ,
        87.        ,  89.        ,  91.15805094,  92.89257099,
        95.39239802,  96.50650778,  99.98896804,  99.62335551,
       105.34812902, 101.32449627, 114.35633594,  81.63376904])

I am not super experienced with digital filters, so maybe I am making a simple error. Why does the python output appear to oscillate at the end of the array? Ultimately, how can I get the python to match the Matlab?

I tried to change the order of the filter in Python with the 'n=XX' parameter in sg.decimate. I see that it seems to behave better with odd values, but it still seems to oscillate at the end of the output array.

I am expecting the Python output to accurately match the Matlab output.

Solution

I ended up following dankal444's advice. The following worked for me:

def FIRdecimate(inputArray, r, nfilt=30):
    if int(r)==1:
        return inputArray.copy()
    b=sg.firwin(nfilt+1, 1/r, window='hamming', pass_zero='lowpass')

    #lead-in points
    inputArrayPadding = 2*inputArray[0] - inputArray[range(nfilt+1, 0, -1)]
    _, zf = sg.lfilter(b, 1, inputArrayPadding, zi=np.zeros(nfilt, dtype=int))
    outputArray, zf = sg.lfilter(b,1,inputArray, zi=zf)

    #lead-out points
    inputArrayPadding=2*inputArray[len(inputArray)-1]-inputArray[range(len(inputArray)-2, len(inputArray)-2*nfilt-4, -1)]
    outputArrayPadding,_ = sg.lfilter(b, 1, inputArrayPadding, zi=zf)
    
    #get start point
    _, gd = sg.group_delay((b,1),8)
       
    #take every rth point
    ndxs = []
    for i in range(int(gd[0]+ 0.75), len(inputArray), r):
        ndxs.append(i)
    
    outputArray = outputArray[ndxs]
    startpoint = r - (len(inputArray)-ndxs[-1])
    outputArray = np.concatenate((outputArray, outputArrayPadding[range(startpoint, startpoint+int(np.ceil(len(inputArray)/r)-len(outputArray))*r-1, r)]))
    return outputArray