0

我正在处理每天以 32hz 出现的大量数据。我们希望将数据过滤到 .5hz（编辑：问题最初指定为 1hz，根据建议进行了更改）并将其采样到 1hz。我正在使用 signal.resample 进行下采样，并使用带有 signal.butter 过滤器的 signal.filtfilt。然而，在执行这两项操作后，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 尖峰由于其他频率泄漏）：
原始数据（4hz 尖峰由于其他频率泄漏）

过滤后：
过滤后

重采样后：
重采样后

编辑：调整到 0.5hz 的过滤器后，会出现同样的问题。我现在想知道我如何显示我的 FFT 是否有问题。我已经包含了我用来显示图表的快速绘图。

4

0 回答 0