параметры фильтра низких частот для пихты - PullRequest
Новая
       17

параметры фильтра низких частот для пихты

6 голосов
/ 11 ноября 2010

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

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

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

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

В другой программе запуск фильтра нижних частот с помощью команд графического интерфейса пользователя приводит к выводу, который имеет аналогичные средние значения для каждого 10-секундного сегмента (10000 точек данных), но который имеет значительно более низкие стандартные отклонения, так что мы практически теряем шум в этом конкретном файле данных и замените его чем-то, что сохраняет среднее значение, показывая более долгосрочные тренды, которые не загрязнены высокочастотным шумом. Диалоговое окно параметров другого программного обеспечения содержит флажок, который позволяет вам выбрать количество коэффициентов, чтобы оно «оптимизировалось на основе размера выборки и частоты выборки». (У меня 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

Ответы [ 2 ]

26 голосов
/ 18 ноября 2010

Отсечка нормализуется до частоты Найквиста, которая составляет половину частоты дискретизации. Таким образом, при 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()

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

1 голос
/ 11 ноября 2010

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

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

...