Implementação de PyWavelets CWT

Question

Eu procuro entender a implementação de PyWavelets da Transformada Wavelet Contínua e como ela se compara à versão mais 'básica' que codifiquei e forneci aqui. Em particular:

1. Como a wavelet integrada, depois a pós-convolução diferenciada, é equivalente à convolução sem int + diff?
2. Como reamostrar a wavelet integrada em resolução crescente é equivalente a alterar sua escala?

---

Implementação ingênua :

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 (preso às partes relevantes) :

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

---

Comparação em $f=1, 4$sinusóides :

Answer 1

Primeiro, consulte "Análise ingênua" na seção abaixo. Em PyWavelets: o algoritmo foi encontrado no Github para derivar de uma implementação MATLAB antiga , mas não fornece detalhes sobre a codificação do kernel wavelet. Então comecei a explorar; vários visuais e anotações de código aqui . Seguem as principais conclusões. - Todos os códigos

Consulte "Análise do PyWavelets" abaixo. Ele está implementando o seguinte - a saber, Eq 4:

---

(1) Como a wavelet integrada, e depois a pós-convolução diferenciada, é equivalente à convolução sem int + diff?

A formulação está na direita, derivando da esquerda, equações; Em vez de integrar continuamente tudo$\mathcal{R}$, a integral é dividida em segmentos que somam $k$. A entrada é considerada constante ao longo$[k, k+1]$ (embora não devesse ser $)$?), pois é discreto, enquanto a wavelet é mantida contínua. Isso permite mover a entrada para fora do integral; visualmente:

Então, como codificamos isso? Está claro pela Eq 4 que precisamos de um difflugar. pywtfaz np.diff(conv), que é basicamente np.diff(coef); por que parece tirar a diferença do coeficiente suposto em Naive?

Na Eq 4, temos a wavelet integrada up-to k + 1, menos up-to k, multiplicada pela wavelet. Isso é equivalente ao produto com a wavelet em b - 1e b, respectivamente (conecte t=k+1e t=k, compare com conectar b=b-1e em vez b=bdisso).

A resposta está então em conv; suponha que ele já tenha sido cortado de forma que seja apenas quando a wavelet estiver pelo menos na metade do sinal. Então, cada ponto em convé o produto da wavelet com uma parte da entrada, em uma diferente tau(ou melhor b). Portanto, dois pontos adjacentes em convsão produtos com wavelets em be b + 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} $$

No entanto, buscamos menos entre t=k+1e t, ou um maior te mais baixo t, que é o negativo de acima, explicando assim a etapa 7 em PW Breakdown.

---

(2) Como reamostrar a wavelet integrada em resolução crescente é equivalente a mudar sua escala?

Veja a figura em (5) em Repartição do PW. pywtparece "escalar" por meio do número de amostras que definem a wavelet, em vez de recomputar a wavelet com um argumento dilatado. ... apesar de, visualmente, este reamostragem claramente faz dilatar a wavelet, esticá-lo em relação à entrada.

Portanto, é de fato uma forma de dimensionamento; mas como ele se compara ao ingênuo? O ingênuo é recalculado de forma mais intuitiva por equações, usando comprimento de wavelet fixo e passando um argumento em escala. Então, o que é melhor? A pergunta oposta é mais clara de responder; considere o que acontece em escalas superiores e inferiores:

1. Ingênuo superior : a wavelet se dilata muito, saindo do quadro de amostragem muito antes de sua cauda zero.
2. pywt higher : wavelet is just fine; o problema está com conv. O comprimento da wavelet é fixado em 1024, portanto, se a entrada for mais curta, as wavelets de escala mais alta nunca podem multiplicar totalmente o sinal. Quanto maior a disparidade, mais a wavelet é "vista" semelhante a "Ingênua mais alto" pelo sinal; isso pode ser visto nos mapas de calor da pergunta diferindo por deslocamentos verticais.
   • A escala máxima da wavelet também é ilimitada; len(j)cresce proporcionalmente com a escala (mas há um problema de normalização, veja aqui ).
3. Ingênuo inferior : tudo bem.
4. pywt lower : wavelet tem resolução pequena.
5. Comprimento de janela fixo vs variável : escalas mais altas devem corresponder a maior dilatação, ou ser diferente de zero, da wavelet em relação à entrada . Ambos conseguem isso, exceto Ingênuo descontinuado para zero abruptamente em escalas mais altas.

Os problemas de escala superior em ambos são resolvidos limitando a escala máxima; se ambos estiverem em suas zonas completamente seguras, a implementação Naive parece preferível devido à resolução mais alta. No entanto, esta "zona segura" para Ingênuo é proibitivamente pequena, se formos considerar "segura" como não saltar para zero na cauda.

Ao contrário, pywto problema de é resolvido mais facilmente aumentando o comprimento máximo da wavelet, o que foi sugerido para livrar-se dos artefatos semelhantes a zíper (que podem ser vistos nos gráficos da pergunta). Isso, no entanto, pioraria o desempenho em maior escala para entradas curtas - mas na prática é improvável que seja menor do que 1024.

---

Análise de PyWavelets :

1. Wavelet, antes da integração, combina exatamente com o blob de código mostrado, que é uma aproximação do Morlet real completo (usado por Naive) assumindo $\sigma > 5$no Wiki .
2. pywtintegra Morlet real via np.cumsum(psi) * step, contabilizando o tamanho do passo diferencial
3. A wavelet integrada,, int_psié reutilizada para todosscales
4. Para cada um scale, o mesmo int_psié reamostrado em resolução crescente, com jum número inteiro variando linearmente entre o mesmo mínimo e máximo (quase) para todos scales:

1. A wavelet resultante para aumentar scale, lado a lado com ingênuo (para mais tarde):

1. O resultado da convolução primeiro diff'd,
2. negado,
3. multiplicado por sqrt(scale),
4. então aparado em ambos os lados, quer (a) incluem apenas partes onde a wavelet é, pelo menos, a meio caminho "dentro do sinal", ou (b) obrigar o seu comprimento a len(x).

---

Repartição ingênua : se não estiver familiarizado com CWT / STFT, recomendo fortemente as partes 1-3 aqui .

Precisamos de (1) tipo de wavelet; (2) comprimento de wavelet; (3) incremento / sobreposição de wavelet. Usaremos Morlet real , 96 amostras e incremento de 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

Visualizado por scale=2:

A diferença entre esta implementação e a implementação ingênua na questão é que o último usa win_inc=1, equivalente a convolução com uma wavelet invertida e aparar o resultado para incluir apenas as partes onde a wavelet estava totalmente "dentro do sinal". Além disso, aqui temos tau("incremento") no loop externo, há escala, mas os dois são equivalentes.

Por último, as formas ao longo da atenuação do timeshift diferem; a implementação da questão apara de forma que apenas as partes da convolução em que a wavelet está pelo menos na metade "dentro da entrada" sejam incluídas; este corta ainda mais onde a wavelet completa está incluída.

---

NÃO SOLVIDO :

1. Por que * sqrt(scale)(etapa 8 pywt)? Eq 4 mostra 1 / sqrt(scale). Eu considerei * scalepor regra de cadeia, exceto que integramos a wavelet numericamente ... Faz sentido ter um tratamento por escala; podemos nos safar sem integrar uma wavelet por escala porque a wavelet por escala é ... a mesma wavelet reamostrada? Resolvido .