Параметры фильтра нижнего пропускания с использованием scipy

Question

Я пытаюсь написать простой фильтр нижних частот, используя scipy, но мне нужна помощь в определении параметров.

У меня есть 3,5 миллиона записей в данных временных рядов, которые необходимо фильтровать, а данные отбираются на 1000 Гц.

Я использую signal.firwin и signal.lfilter из библиотеки scipy.

Параметры, которые я выбираю в приведенном ниже коде, не фильтруют мои данные вообще. Вместо этого приведенный ниже код просто создает то, что графически выглядит как одни и те же точные данные, за исключением искажений фазы времени, которые сдвигают график вправо чуть менее чем на 1000 точек данных (1 секунду).

В другом программном обеспечении, запускающем фильтр нижних частот с помощью графических команд пользовательского интерфейса, выводится вывод, который имеет аналогичные средства для каждого сегмента 10 секунд (10 000 точек данных), но имеет существенно более низкие стандартные отклонения, так что мы по существу теряем шум в этом конкретном файле данных и заменить его чем-то, что сохраняет среднее значение, показывая долгосрочные тенденции, которые не загрязняются более высоким уровнем шума. Диалоговое окно параметров другого программного обеспечения содержит флажок, который позволяет вам выбирать количество коэффициентов, чтобы оно «оптимизировалось на основе размера выборки и частоты дискретизации». (Mine 3,5 миллиона образцов, собранных в 1000 герц, но я хотел бы функцию, которая использует эти входы в качестве переменных.)

* Может кто-нибудь показать мне, как настроить этот код таким образом, что она удаляет все частоты выше 0,05 хз? * Я хотел бы видеть гладкие волны на графике, а не только временные искажения того же идентичного графа, что и сейчас, из кода ниже.

class FilterTheZ0(): 
    def __init__(self,ZSmoothedPylab): 
     #------------------------------------------------------ 
     # Set the order and cutoff of the filter 
     #------------------------------------------------------ 
     self.n = 1000 
     self.ZSmoothedPylab=ZSmoothedPylab 
     self.l = len(ZSmoothedPylab) 
     self.x = arange(0,self.l) 
     self.cutoffFreq = 0.05 

     #------------------------------------------------------ 
     # Run the filter 
     #------------------------------------------------------ 
     self.RunLowPassFIR_Filter(self.ZSmoothedPylab, self.n, self.l 
             , self.x, self.cutoffFreq) 

    def RunLowPassFIR_Filter(self,data, order, l, x, cutoffFreq): 
     #------------------------------------------------------ 
     # Set a to be the denominator coefficient vector 
     #------------------------------------------------------ 
     a = 1 
     #---------------------------------------------------- 
     # Create the low pass FIR filter 
     #---------------------------------------------------- 
     b = signal.firwin(self.n, cutoff = self.cutoffFreq, window = "hamming") 

     #--------------------------------------------------- 
     # Run the same data set through each of the various 
     # filters that were created above. 
     #--------------------------------------------------- 
     response = signal.lfilter(b,a,data) 
     responsePylab=p.array(response) 

     #-------------------------------------------------- 
     # Plot the input and the various outputs that are 
     # produced by running each of the various filters 
     # on the same inputs. 
     #-------------------------------------------------- 

     plot(x[10000:20000],data[10000:20000]) 
     plot(x[10000:20000],responsePylab[10000:20000]) 
     show() 

     return 

Answer 1

Единицы частоты отсечки, вероятно, [0,1], где 1.0 эквивалентно FS (частота дискретизации). Поэтому, если вы действительно имеете в виду 0,05 Гц и FS = 1000 Гц, вы должны пройти cutoffFreq/1000. Вам может понадобиться более длинный фильтр, чтобы получить такую ​​низкую отсечку.

(BTW вы передаете какие-то аргументы, но затем с помощью атрибутов объекта вместо этого, но я не вижу, что введение каких-либо очевидных ошибок все же ...)

Answer 2

Граничная нормирована на частоте Найквиста, которая составляет половину частота выборки. Таким образом, при FS = 1000 и FC = 0,05 вы хотите, чтобы обрезание = 0,05/500 = 1e-4.

from scipy import signal 

FS = 1000.0           # sampling rate 
FC = 0.05/(0.5*FS)         # cutoff frequency at 0.05 Hz 
N = 1001            # number of filter taps 
a = 1            # filter denominator 
b = signal.firwin(N, cutoff=FC, window='hamming') # filter numerator 

M = FS*60           # number of samples (60 seconds) 
n = arange(M)          # time index 
x1 = cos(2*pi*n*0.025/FS)       # signal at 0.025 Hz 
x = x1 + 2*rand(M)         # signal + noise 
y = signal.lfilter(b, a, x)       # filtered output 

plot(n/FS, x); plot(n/FS, y, 'r')     # output in red 
grid() 

Выходной фильтр задерживается полсекунды (фильтр по центру на кране 500). Обратите внимание, что смещение постоянного тока, добавленное шумом, сохраняется фильтром нижних частот. Кроме того, 0,025 Гц находится в пределах диапазона прохождения, поэтому выходное колебание от пика до пика составляет приблизительно 2.