matzner_bar-gad_PLoS_2015.py 文件源码

def bandfilter(data,fa=None,fb=None,Fs=1000.,order=4,zerophase=True,bandstop=False):
        N = len(data)
        assert len(shape(data))==1
        padded = zeros(2*N,dtype=data.dtype)
        padded[N/2:N/2+N]=data
        padded[:N/2]=data[N/2:0:-1]
        padded[N/2+N:]=data[-1:N/2-1:-1]
        if not fa==None and not fb==None:
            if bandstop:
                b,a  = butter(order,array([fa,fb])/(0.5*Fs),btype='bandstop')
            else:
                b,a  = butter(order,array([fa,fb])/(0.5*Fs),btype='bandpass')
        elif not fa==None:
            # high pass
            b,a  = butter(order,fa/(0.5*Fs),btype='high')
            assert not bandstop
        elif not fb==None:
            # low pass
            b,a  = butter(order,fb/(0.5*Fs),btype='low')
            assert not bandstop
        else:
            assert 0
        if zerophase:
            return filtfilt(b,a,padded)[N/2:N/2+N]
        else:
            return lfilter(b,a,padded)[N/2:N/2+N]
        assert 0

    # now try it with a broad spectrum