6

scipy を使用して単純なローパス フィルターを作成しようとしていますが、パラメーターの定義に助けが必要です。

フィルター処理が必要な時系列データに 350 万件のレコードがあり、データは 1000 Hz でサンプリングされています。

scipy ライブラリの signal.firwin と signal.lfilter を使用しています。

以下のコードで選択しているパラメーターは、データをまったくフィルタリングしません。代わりに、以下のコードは、グラフを 1000 データ ポイント (1 秒) 弱だけ右にシフトする時間位相の歪みを除いて、グラフィック的にはまったく同じデータのように見えるものを単純に生成します。

別のソフトウェア プログラムでは、グラフィカル ユーザー インターフェイス コマンドを介してローパス モミ フィルターを実行すると、10 秒 (10,000 データ ポイント) セグメントごとに同様の平均値を持つ出力が生成されますが、標準偏差が大幅に低下するため、この特定のノイズが本質的に失われます。より高い周波数のノイズによって汚染されていない長期的な傾向を示しながら、平均値を保持するものに置き換えます。他のソフトウェアのパラメーター ダイアログ ボックスには、「サンプル サイズとサンプリング周波数に基づいて最適化する」ように係数の数を選択できるチェック ボックスが含まれています。(私のものは 1000 Hz で収集された 350 万のサンプルですが、これらの入力を変数として使用する関数が必要です。)

* 0.05 Hz を超えるすべての周波数を削除するように、以下のコードを調整する方法を誰か教えてもらえますか? * 以下のコードから取得している同じ同一のグラフの時間の歪みだけでなく、グラフで滑らかな波を見たいです。

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
4

2 に答える 2

26

カットオフは、サンプリング レートの半分であるナイキスト周波数に正規化されます。したがって、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()

出力 フィルター出力は 0.5 秒遅延します (フィルターはタップ 500 を中心にしています)。ノイズによって追加された DC オフセットは、ローパス フィルターによって保持されることに注意してください。また、0.025 Hz は十分に合格範囲内にあるため、ピーク間の出力スイングは約 2 です。

于 2010-11-18T05:16:15.990 に答える
1

カットオフ周波数の単位はおそらく [0,1) で、1.0 は FS (サンプリング周波数) に相当します。したがって、本当に 0.05Hz と FS=1000Hz を意味する場合は、パスする必要がありますcutoffFreq / 1000。このような低いカットオフを得るには、より長いフィルターが必要になる場合があります。

（ところで、いくつかの引数を渡していますが、代わりにオブジェクト属性を使用していますが、明らかなバグがまだ導入されているとは思いません...）

于 2010-11-11T08:38:20.587 に答える