Python中对同频数据的重采样和滤波

Question

我在处理每天以32hz频率输入的大量数据。我们希望将数据过滤到0.5hz（编辑：问题最初指定为1hz，根据建议更改），并将其向下采样到1hz。我正在使用信号重采样用于下采样和信号过滤器用一个信号。黄油过滤器。然而，在执行这两项之后，FFT显示信号下降约0.16hz。在

你过滤的顺序比下采样重要吗？这和程序有关吗？我不明白这些方法对吗？在

我包括了我认为相关的数据

desiredFreq = 1 * 60 *60 *26

# untouched data
quickPlots(x,32) # 32 hz

# filtered
x = ft.butterLowPass(x, .5, 32)
quickPlots(x,32) # 32 hz

# resampled
x = ft.resampleData(x, desiredFreq)
quickPlots(x,1) # should be 1 hz


def butterLowPass(data,desired,original,order = 5):
    """ Creates the coefficients for a low pass butterworth
    filter. That can change data from its original sample rate
    to a desired sample rate

Args: 
----
    data (np.array): data to be filtered
    desirerd (int): the cutoff point for data in hz
    original (int): the initial sampling rate in hz

Returns:
-------
    np.array: of the data after it has been filtered

Note:
----
    I find that the filter will make the data all 0's if the order
is too high. Not sure if this is my end or scipy.signal. So work with
CAUTION

References:
----------
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.butter.html

https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.filtfilt.html

https://en.wikipedia.org/wiki/Butterworth_filter

https://en.wikipedia.org/wiki/Low-pass_filter

https://en.wikipedia.org/wiki/Nyquist_rate
"""
nyq = .5 * original
cutoff = desired / nyq
b, a = signal.butter(order, cutoff, btype = 'lowpass', analog = False)

return signal.filtfilt(b, a, data, padlen =10)

def resampleData(data, desired):
    """ Takes in a set of data and resamples the data at the 
    desired frequency.



Args:
----
    data (np.array): data to be operated on
    desired (int): the desired sampled points 
        over the dataset

Returns:
-------
    np.array: of the resamples data

Notes:
-----
    Since signal.resample works MUCH faster when the len of 
data is a power of 2, we pad the data at the end with null values
and then discard them after.
"""
nOld = len(data)    
thispow2 = np.floor(np.log2(nOld))
nextpow2 = np.power(2, thispow2+1)

data = np.lib.pad(
        data, 
        (0,int(nextpow2-nOld)), 
        mode='constant',
        constant_values = (0,0)
        )

nNew = len(data)
resample = int(np.ceil(desired*nNew/nOld)) 


data = signal.resample(data, resample)    
data = data[:desired]
return data

def quickPlots(data,fs):
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from scipy import signal
    import time
    fft,pxx=signal.welch(data,fs)
    plt.plot(data)
    plt.show()
    plt.loglog(fft,pxx)
    plt.show()

FFT图片：

原始数据（由于其他频率泄漏，4hz峰值）：

filterfilt之后：

重采样后：

编辑：调整到0.5hz的滤波器后，同样的问题出现了。我现在想知道是不是我显示FFT的方式有问题。我包括了我用来显示图形的快速绘图。在