PyWavelets CWT 구현

Sep 28 2020

저는 PyWavelets 의 Continuous Wavelet Transform 구현 을 이해 하고 여기에서 코딩하고 제공 한보다 ' 기본'버전과 비교하는 방법 을 이해하려고합니다 . 특히:

통합 된 웨이블릿은 int + diff없이 convolving하는 것과 같은 방식으로 post-convolving 차이가 나나요?
스케일을 변경하는 것과 동일한 해상도로 통합 웨이블릿을 리샘플링하는 방법은 무엇입니까?

순진한 구현 :

def cwt(x, win_len=None, win='morlet', real=False):
    N = len(x)
    win_len = win_len or N // 8

    scales = _scales(N, nv=32)
    coef = np.zeros((len(scales), N), dtype='complex128')
    kernel = morlet_kernel(win_len)
    wl2 = win_len // 2

    for i, scale in enumerate(scales):
        coef[i, :] = np.convolve(x, kernel(scale)[::-1])[wl2:-(wl2 - 1)]
    return coef

PyWavelets (관련 부분으로 잘림) :

def pywt_cwt(data, scales):
    out = np.empty((np.size(scales),) + data.shape)
    int_psi, x = integrate_wavelet(wavelet='morl', precision=10)
    
    for i, scale in enumerate(scales):
        step = x[1] - x[0]
        j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step)
        j = j.astype(int)  # floor
        if j[-1] >= int_psi.size:
            j = np.extract(j < int_psi.size, j)
        int_psi_scale = int_psi[j][::-1]

        conv = np.convolve(data, int_psi_scale)
        coef = - np.sqrt(scale) * np.diff(conv, axis=-1)
        
        d = (coef.shape[-1] - data.shape[-1]) / 2.
        coef = coef[..., floor(d):-ceil(d)]
        out[i, ...] = coef
    return out

비교 $f=1, 4$정현파 :

답변

OverLordGoldDragon Sep 28 2020 at 20:30

먼저 아래 섹션의 "순진한 분석"을 참조하십시오. PyWavelets에 대해 : 알고리즘은 Github 에서 이전 MATLAB 구현 에서 유래 한 것으로 나타 났지만 웨이블릿 커널 코딩에 대한 세부 정보는 제공하지 않습니다. 그래서 저는 탐험을 시작했습니다. 여기에 다양한 비주얼과 코드 주석이 있습니다 . 주요 결과는 다음과 같습니다. - 모든 코드

아래의 "PyWavelets 분석"을 참조하십시오. 다음을 구현하고 있습니다. 즉, Eq 4입니다.

(1) int + diff없이 convolving과 동등하게 통합 된 웨이블릿과 차이가있는 post-convolving은 어떻게됩니까?

공식은 오른쪽에 있으며 왼쪽에서 파생 된 방정식입니다. 모든 것을 지속적으로 통합하는 대신$\mathcal{R}$, 적분은 다음과 같은 세그먼트로 분할됩니다. $k$. 입력이 일정하다고 가정합니다.$[k, k+1]$ (그럴 수는 없지만 $)$?), 이산 적이므로 잔물결은 연속적으로 유지됩니다. 이는 입력을 적분 외부로 이동할 수있게합니다. 시각적으로 :

그럼 어떻게 코딩할까요? 식 4에서 우리에게 diff어딘가 가 필요하다는 것이 분명 합니다. pywt하지 np.diff(conv), 그것은 기본적으로 np.diff(coef); Naive에서 예상 계수의 차이가 나타나는 이유는 무엇입니까?

식 4에서 우리는 통합 웨이블릿 up-to k + 1, minus up-to k, 웨이블릿을 곱했습니다. 이는 웨이블릿이 각각 및에 있는 제품과 동일 합니다 ( 및 을 연결 하고 대신 연결하는 것과 비교).b - 1bt=k+1t=kb=b-1b=b

대답은 다음 안에 있습니다 conv. 웨이블릿이 적어도 신호 내부의 중간에있을 때만 해당되도록 이미 트리밍되었다고 가정합니다. 그런 다음 모든 지점은 conv다른 tau(또는 오히려 b) 에서 입력의 일부와 웨이블릿의 곱입니다 . 따라서의 인접한 두 지점은 및에 conv웨이블릿이있는 제품입니다 .bb + 1

$$ \begin{align} \text{conv} &= [s \cdot \psi_{\text{int}}(t - 0),\ s \cdot \psi_{\text{int}}(t - 1),\ ...] \\ \Rightarrow \text{DIFF}(\text{conv}) &= [s \cdot (\psi_{\text{int}}(t - 1) - \psi_{\text{int}}(t - 0)),\ ...] \end{align} $$

그러나 사이 마이너스 추구 t=k+1하고 t, 또는 더 t낮출 t따라서 PW 고장 7 단계를 설명하기 위의 네가티브이다.

(2) 스케일을 변경하는 것과 동일한 해상도로 통합 웨이블릿을 리샘플링하는 방법은 무엇입니까?

PW 분석의 (5) 아래 그림을 참조하십시오. pywt확장 된 인수로 웨이블릿을 재 계산하는 대신 웨이블릿을 정의하는 샘플 수를 통해 "규모"하는 것처럼 보입니다. ... 비록 시각이 명확 리샘플링하여 수행 입력하는 것이 상대적으로 스트레칭 웨이블릿 팽창.

그래서 그것은 실제로 확장의 한 형태입니다. 그러나 그것은 순진한 것과 어떻게 비교됩니까? Naive는 고정 된 웨이블릿 길이를 사용하고 스케일링 된 인수를 전달하여 방정식별로 더 직관적으로 다시 계산됩니다. 그렇다면 어느 것이 더 낫습니까? 정반대의 질문이 더 명확합니다. 더 높고 낮은 규모에서 일어나는 일을 고려하십시오.

더 높은 순진함 : 잔물결이 너무 많이 확장되어 제로 테일보다 훨씬 앞서 샘플링 프레임을 빠져 나갑니다.
pywt 더 높음 : 잔물결은 괜찮습니다. 문제는 conv. 웨이브 렛 길이는 1024로 고정되어 있으므로 입력이 더 짧으면 더 높은 스케일의 웨이브 렛이 신호를 완전히 곱할 수 없습니다. 불일치가 클수록 신호에 의해 "Naive high"와 유사한 웨이블릿이 더 많이 "보여집니다". 이것은 수직 이동에 따라 다른 질문의 히트 맵에서 볼 수 있습니다.
- 잔물결의 최대 규모도 무제한입니다. len(j)규모에 비례하여 증가합니다 (하지만 정규화 문제가 있습니다. 여기 참조 ).
Naive Lower : 모두 좋습니다.
pywt lower : 잔물결의 해상도가 작습니다.
고정 대 가변 창 길이 : 더 높은 스케일은 입력에 비해 웨이블릿의 더 큰 확장 또는 0이 아닌 것에 해당해야합니다 . Naive가 더 높은 스케일에서 갑자기 제로화를 중단하는 것을 제외하고는 둘 다이를 수행합니다.

두 가지 모두에서 더 높은 척도 문제는 최대 척도를 제한하여 해결됩니다. 둘 다 완전히 안전한 영역에 있다면 Naive 구현이 더 높은 해상도로 인해 선호되는 것 같습니다. 그러나 "안전"을 꼬리에서 0으로 점프하지 않는 것으로 간주한다면 Naive의 "안전 영역"은 엄청나게 작습니다.

반대로, pywt의 문제는 지퍼와 같은 인공물을 제거하기 위해 제안 된 최대 웨이블릿 길이를 늘림으로써 더 쉽게 해결됩니다 (질문의 플롯에서 볼 수 있음). 그러나 이는 짧은 입력에 대해 더 높은 규모의 성능을 악화 시키지만 실제로는 1024보다 짧을 것 같지 않습니다.

PyWavelets 분석 :

통합 전에 Wavelet은 표시된 코드 blob과 정확히 일치하며, 이는 전체 실제 Morlet (Naive에서 사용)의 근사치입니다. $\sigma > 5$에서 위키 .
pywt를 통해 실제 Morlet을 통합 np.cumsum(psi) * step하여 차동 단계 크기를 설명합니다.
통합 웨이블릿 int_psi은 모두에게 재사용됩니다.scales
각각 scale에 대해 동일한 값 int_psi이 증가하는 해상도 로 리샘플링 되며 j정수는 동일한 최소값과 최대 값 (거의) 사이에서 선형으로 변합니다 scales.

순진함 scale과 나란히 증가하는 웨이블릿 (나중을 위해) :

컨볼 루션 결과 첫 번째 diff'd,
부정,
곱하기 sqrt(scale),
그런 다음 (a) 웨이블릿이 "신호 내부"의 절반 이상인 부분 만 포함하거나 (b) 길이를 len(x).

Naive Breakdown : CWT / STFT에 익숙하지 않은 경우 여기에서 1-3 부를 강력히 권장 합니다 .

우리는 (1) 잔물결 유형이 필요합니다. (2) 잔물결 길이; (3) 웨이블릿 증가 / 중첩. 실제 Morlet , 96 개 샘플을 사용하고 20 씩 증가합니다.

win_len = 96  # "win" == window == wavelet
win_inc = 20
n_wins = (len(x) - win_len) // win_inc + 1

coef = np.zeros((n_wins, len(scales)), dtype='complex128')
kernel = morlet_kernel(win_len)

for tau in range(n_wins):
    start = tau * win_inc
    end   = start + win_len
    coef[tau, :] = _transform(x[start:end], kernel, scales)

def _transform(x, kernel, scales):
    coef = np.zeros(len(scales), dtype='complex128')
    for i, scale in enumerate(scales):
        psi = np.conj(kernel(scale))
        coef[i] = np.sum(x * psi / np.sqrt(scale))
    return coef

시각화 scale=2:

질문에서이 구현과 순진한 구현의 차이점은 후자는 win_inc=1뒤집힌 웨이블릿으로 컨 볼빙하고 웨이블릿이 완전히 "신호 내부"인 부분 만 포함하도록 결과를 트리밍하는 것과 동일합니다. 또한 여기 tau바깥 쪽 루프에 ( "증가") 스케일이 있지만 둘은 동일합니다.

마지막으로, 타임 시프트 희미한 모양이 다릅니다. 질문의 구현은 웨이블릿이 적어도 "입력 내부"의 중간에있는 컨볼 루션의 일부만 포함되도록 트리밍됩니다. 이것은 전체 잔물결이 포함 된 곳으로 더 잘 립니다.