2013-04-30 6 views
0

저는 파이썬에서 컷오프 주파수가 0.7-4Hz 인 128 포인트 해밍 윈도우로 밴드 패스 필터를 얻으려고합니다. 이미지에서 내 신호에 대한 샘플을 얻습니다. (1 샘플 = 1 이미지). fps는 종종 변경됩니다.파이썬의 밴드 패스 필터

어떻게 이것을 파이썬으로 할 수 있습니까? 나는 이것을 읽었다 : http://mpastell.com/2010/01/18/fir-with-scipy/ 그러나 나는 혼동하는 firwin을 발견한다. 이 변수 fps로 어떻게 이것을 할 수 있습니까?

+0

샘플링 속도는 어떻게됩니까? 알지 못하면 디지털 컷오프 주파수를 지정할 수 없습니다. – sizzzzlerz

+0

글쎄 초당 프레임 수와 같으니 까. 그것은 끊임없이 변화합니다 ...20fps ~ 60fps 범위 ... 샘플링 주파수가 변경됩니다. – Ojtwist

+0

필터링하려는 데이터가 일정한 속도로 샘플링되지 않는다는 말입니까? 이 경우 FIR 필터 나 다른 디지털 필터를 사용할 수 없습니다. 일정하고 일관된 샘플링 속도가 필요합니다. – sizzzzlerz

답변

2

일관성없는 샘플링 속도로 데이터를 필터링하려고 시도하는 것은 매우 어렵습니다 (불가능합니까?). 그래서 당신이하고자하는 것은 다음과 같습니다.

  1. 고정 된 샘플 속도로 새로운 신호를 생성하십시오. 고정 샘플 속도는 최대 샘플 속도 이상이어야합니다. 새 샘플을 가져올 위치를 나타내는 새로운 "그리드"를 설정하고 기존 데이터의 값을 보간하여이 작업을 수행하십시오. 얼마나 정확하게해야하는지에 따라 다양한 보간 방법이 있습니다. 선형 보간은 나쁜 출발점이 아니지만 실제로하는 일에 달려 있습니다. 확실하지 않으면 https://dsp.stackexchange.com/에 문의하십시오.

  2. 일단 작업을 완료하면 링크 된 게시물에 설명 된 것과 같이 샘플이 고르게 배치되므로 표준 신호 처리 방법을 신호에 적용 할 수 있습니다.

  3. 필요한 경우 원래 샘플 위치를 다시 얻으려면 다시 보간해야 할 수도 있습니다.

데이터를 분석하려는 경우 Lomb Periodigram에 관심이있을 수 있습니다. 데이터를 밴드 패스 한 다음 분석하면, Lomb Periodigram을 사용하고 관련 주파수 만 보거나 원하는 결과를 가중시킵니다. (숫자 레시피 시리즈를 참조하십시오.) 13.8 절은 "wikipedia 페이지보다 우호적 인 소개 인 것처럼 보이는 편평하지 않은 데이터의 스펙트럼 분석"이라고합니다.

+0

이것은 들리는데, 이것에 대해 살펴 보겠습니다. – Ojtwist

+0

FYI : SciPy에는 Lomb-Scargle의 주기도 인 scipy.signal.lombscargle이 구현되어 있습니다. –

8

scipy.signal.firwin 또는 scipy.signal.firwin2 함수를 사용하여 밴드 패스 FIR을 만들 수 있습니다 필터. scipy.signal.remez

다음 코드는 밴드 패스 FIR 필터를 만들기위한 몇 가지 편리한 래퍼를 제공합니다. 이것들을 사용하여 질문에서 요구 한 숫자에 해당하는 대역 통과 필터를 만듭니다. 이것은 샘플링이 균일하게 수행된다고 가정합니다. 샘플링이 일정하지 않으면 FIR 필터가 적합하지 않습니다.

from scipy.signal import firwin, remez, kaiser_atten, kaiser_beta 

# Several flavors of bandpass FIR filters. 

def bandpass_firwin(ntaps, lowcut, highcut, fs, window='hamming'): 
    nyq = 0.5 * fs 
    taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False, 
        window=window, scale=False) 
    return taps 

def bandpass_kaiser(ntaps, lowcut, highcut, fs, width): 
    nyq = 0.5 * fs 
    atten = kaiser_atten(ntaps, width/nyq) 
    beta = kaiser_beta(atten) 
    taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False, 
        window=('kaiser', beta), scale=False) 
    return taps 

def bandpass_remez(ntaps, lowcut, highcut, fs, width): 
    delta = 0.5 * width 
    edges = [0, lowcut - delta, lowcut + delta, 
      highcut - delta, highcut + delta, 0.5*fs] 
    taps = remez(ntaps, edges, [0, 1, 0], Hz=fs) 
    return taps 


if __name__ == "__main__": 
    import numpy as np 
    import matplotlib.pyplot as plt 
    from scipy.signal import freqz 

    # Sample rate and desired cutoff frequencies (in Hz). 
    fs = 63.0 
    lowcut = 0.7 
    highcut = 4.0 

    ntaps = 128 
    taps_hamming = bandpass_firwin(ntaps, lowcut, highcut, fs=fs) 
    taps_kaiser16 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.6) 
    taps_kaiser10 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.0) 
    remez_width = 1.0 
    taps_remez = bandpass_remez(ntaps, lowcut, highcut, fs=fs, 
           width=remez_width) 

    # Plot the frequency responses of the filters. 
    plt.figure(1, figsize=(12, 9)) 
    plt.clf() 

    # First plot the desired ideal response as a green(ish) rectangle. 
    rect = plt.Rectangle((lowcut, 0), highcut - lowcut, 1.0, 
         facecolor="#60ff60", alpha=0.2) 
    plt.gca().add_patch(rect) 

    # Plot the frequency response of each filter. 
    w, h = freqz(taps_hamming, 1, worN=2000) 
    plt.plot((fs * 0.5/np.pi) * w, abs(h), label="Hamming window") 

    w, h = freqz(taps_kaiser16, 1, worN=2000) 
    plt.plot((fs * 0.5/np.pi) * w, abs(h), label="Kaiser window, width=1.6") 

    w, h = freqz(taps_kaiser10, 1, worN=2000) 
    plt.plot((fs * 0.5/np.pi) * w, abs(h), label="Kaiser window, width=1.0") 

    w, h = freqz(taps_remez, 1, worN=2000) 
    plt.plot((fs * 0.5/np.pi) * w, abs(h), 
      label="Remez algorithm, width=%.1f" % remez_width) 

    plt.xlim(0, 8.0) 
    plt.ylim(0, 1.1) 
    plt.grid(True) 
    plt.legend() 
    plt.xlabel('Frequency (Hz)') 
    plt.ylabel('Gain') 
    plt.title('Frequency response of several FIR filters, %d taps' % ntaps) 

    plt.show() 

다음은 스크립트에 의해 생성 된 플롯입니다. 물론 스크립트를 로컬에서 실행하는 것이 훨씬 더 유용하므로 세부 사항을 확대 할 수 있습니다.

plot of frequency responses

0

또 다른 옵션은 이전에 필터링 ("리딩"와 같은) 일정한 샘플 속도로 데이터를 변환하는 (비동기) 샘플 레이트 변환 될 것이다. 물론 이것은 샘플 속도를 알고있는 경우에만 작동하며 데이터를 필터링해야하고 스펙트럼을 추정 할 필요가없는 경우에만 유용합니다. 이 목적을 위해, 예를 들어,. scipy.interpolate의 InterpolatedUnivariateSpline은 프레임 단위로 적용 할 수 있습니다. 매우 빠릅니다.