PyWavelets CWT: normalização? Vs Scipy?
Relacionado . A equação que está sendo implementada normaliza por sqrt(1 / scale):
$$ C_{a, b} = \frac{1}{\sqrt{a}} \sum_k s(k)\left( \int_{-\infty}^{k+1} \overline{\psi \left(\frac{t-b}{a} \right) }dt - \int_{-\infty}^{k} \overline{\psi \left(\frac{t-b}{a} \right) }dt \right) $$
então por que pywt._cwtfazer sqrt * (escala) ? Além disso, como está normalizando (L1, L2, etc)?
Edit : Depois de revisar o cwt de Scipy , percebi que está implementando exatamente o que descrevi como "ideal" na resposta, portanto, é o mesmo que comparar os dois (exceto alguns detalhes óbvios). Não sei em que cenários o algoritmo de pywt é melhor, mas geralmente - vitória fácil para Scipy . Abri um PR para abordar a falha da chave de pywt.
Edit2 : Encontrou uma borda aparente de pywt sobre scipy no lugar menos esperado: escalas baixas . Endereçado aqui .
Edit 3 : depois de consertar o low precisionin pywt, eu realmente achei pywtsuperior scipyem uma variedade de sinais sintéticos, mas squeezepysempre superei os dois. Testes aqui , com mais vindo.
Respostas
Vamos seguir a matemática desde a incubação até o parto. Ele começa com psium reescalonado morlet2(como mostrado anteriormente ) em uma escala$a=64$, e $\sigma=5$:
$$ \psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2} $$
é integrado e normalizado L1 : - (consulte a ressalva 2 abaixo)
$$ \psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3} $$
A norma está errada por um fator escalar, mas o que importa é o mesmo para diferentes$a$; você pode verificar se é assim (para grandes$a$, aumentar os limites de integração). Observe que é a wavelet que é normalizada por L1, não a wavelet integrada$\psi_\text{int}(t)$; a distinção é essencial, como veremos.
$\psi$ é então convolvido com o sinal $s$, e finito diferenciado e negado; a partir de aqui , isto é:
$$ -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4} $$
Inspecionando de perto, isso é equivalente a desfazer a integração anterior, com certas diferenças numéricas (abordadas posteriormente), então ficamos com
$$ \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5} $$
Por último * sqrt(scale):
$$ \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6} $$
que é, portanto, em última análise, L2-normalizada . Podemos verificar isso em um sinal de teste; L2 é tendencioso na representação de amplitude / potência em direção às frequências mais baixas, suprimindo as frequências altas. Se acima estiver correto, não apenas podemos ver isso, mas também aplicar uma correção removendo * sqrt(scale)do código:
Advertências específicas da implementação discutidas a seguir.
Quão confiável é a normalização entre escalas? - Não é ótimo (- code2 ):
À direita, vemos o melhor caso, usando um número limitado de amostras, à esquerda, a implementação da reamostragem. Claramente, devido à resolução inerentemente pobre de um pequeno número de amostras, nossas normas não funcionam bem para escalas baixas - mas a variante recomputada pelo menos mantém uma norma média aproximadamente constante entre elas.
Mais interessante, o reamostrado de repente parece quadrático após $\log_2(a)=6$; qual é o problema? Lembre-se de que "começa ... em escala$a=64$"; isso é quando reamostramos com mais amostras do que no original ( int_psi):
de quem diffé
- Em 1LR vemos
int_psinascale=256é reamostragem os mesmos valores deint_psimúltiplas vezes, que se manifesta como zeros emdiff(2LR). - Em seguida, elevamos ao quadrado os valores para ter uma ideia da norma L2 e comparamos com os valores normalizados L2 em 3LR; encontramos valores normados soprando valores não normados da água. Isso é excessivo ou apropriado?
- Compare com recalculado na mesma escala; em 4LR, pelo menos vemos a forma recomputada - não é mais desprezível. Isso sugere que a reamostragem está superestimando a amplitude da wavelet e, como no cálculo normal, isso é significativamente compensado pela metade dos valores sendo zero devido a
diff, então a superestimação deve ser grande.
Que ótimo? Podemos comparar as amplitudes das ondas de pico entre as escalas do não normalizado reamostrado diff(int_psi_256)e recomputado w( morlet2/scale); normatizar mudanças tanto pelo mesmo fator ( *sqrt(scale)), não fará diferença.
Com certeza, novamente em $a=64$vemos um comportamento crítico; a amplitude de pico da reamostragem satura no mínimo, enquanto a recomputada continua a diminuir. A diferença em$a=256$é mais do que grande o suficiente (~ 12db = 4x, se ao quadrado -> 16x) para compensar os zeros diff, daí porque a norma L2 cresce.
Mas observe que a norma L1 permanece bem comportada em escalas mais altas. Para descobrir por que, devemos comparar precisamente a proporção de zeros com o fator de superestimação. Já sabemos que o último é 4x; para o primeiro, peguei uma região não plana (caso contrário, são zeros) e encontrei a proporção de zeros para o número de amostras: 0,75 . O que significa que 1 em 4 valores são diferentes de zero e cada um é 4x maior ... a proporção exata necessária para manter a norma L1.
Essa proporção de zeros não é muito clara nos gráficos acima - então, em vez disso, os pontos:
Se diffdesfaz cumsum, por que ambos?
Boa pergunta. Na verdade, ao me livrar de ambos, descobri que o efeito de degradação da qualidade (zíper) desapareceu:
Mas algo mais também desaparece - as escalas de alta frequência (à esquerda) e o sinal vira (desfiz isso removendo menos em coef = -). Eu parcialmente compenso isso substituindo * sqrt(scale)por 1 / sqrt(scale)(à direita), mas isso nos leva ao caso normalizado L2 ( 1 / scaleé um exagero), e é uma modificação sem base à luz de tudo o que foi encontrado até agora. Então, isso pode ser feito? Claro, pode-se implementar sem a int-diffdança, mas não vou explorar isso aqui.
Advertência: norma efetiva
Recall, convé aparado para len(x), de modo que toda a wavelet pode não participar do cálculo, enquanto os cálculos acima assumem o contrário. A suposição é mais verdadeira quanto mais longa x- mas eis o que acontece se len(x) == 512, com a mesma wavelet precision=10(2 ^ 10 = 1024 int_psi), em comparação com a wavelet recomputada:
Uma bagunça! Ambos distorcem significativamente, possivelmente implicando que é inerente a escalas grandes o suficiente fixadas em sinais curtos. Uma solução simples é zerar o sinal, mas eu não atestaria isso sem entender as advertências (o que não faço para wavelets, mas faço para DFT ).
Além disso, observe que a norma não é apenas por escala, mas também por turno ($b$); menos da wavelet está no sinal quando ele está entrando, em oposição a quando seu centro está no centro do sinal.
Mas mantendo-o por escala, confesso que trapacei no exemplo anterior com $f=4$; inicialmente eu escolhi$f=1$, mas o gráfico da norma L1 não parecia bom e eu não entendi por quê. Como agora eu faço, aqui estão L1 e L2 por escala para$f=1,64$:
Compare com o gráfico anterior; observe como os modos da esquerda (escalas baixas) de ambas as normas não são afetados, apenas os da escala alta (baixa frequência) são agora muito maiores (L2) e menores (L1), exatamente como esperado.
Advertência 2: * stepnão faz a normalização de L1
Em vez disso, ele resulta da diffondulação integrada reamostrada; o * stepúnico define a escala inicial correta e pode ser omitido sem perder a norma. Isso não acontece com a recomputação da wavelet ou com a reamostragem da wavelet não integrada.
(Observe que recomputado scipy.signal.morlet2é L2-normalizado, mas apenas porque é explicitamente normalizado; o diff-resample- inté " autonormalizado ").