PyWavelets CWT: normalisation? Vs Scipy?
Connexes . L'équation implémentée se normalise par 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) $$
alors pourquoi ne pywt._cwtpas * sqrt (échelle) ? Aussi, comment se normalise-t-il (L1, L2, etc.)?
Edit : Après avoir examiné le cwt de Scipy , j'ai réalisé qu'il implémentait exactement ce que j'ai décrit comme «idéal» dans la réponse, donc c'est la même chose que de comparer les deux (sauf quelques détails évidents). Je ne sais pas dans quels scénarios l'algorithme de pywt est meilleur, mais en général - victoire facile pour Scipy . J'ai ouvert un PR pour corriger la faille clé de pywt.
Edit2 : Trouvé un bord apparent de pywt sur scipy à l'endroit le moins attendu: les échelles basses . Adressé ici .
Edit 3 : après avoir corrigé le low precisionin pywt, j'ai trouvé en fait pywtsupérieur à scipyune variété de signaux synthétiques, mais je squeezepybat toujours les deux. Tests ici , avec plus à venir.
Réponses
Suivons le calcul de l'incubation à la livraison. Il commence par psi, un redimensionné morlet2(comme indiqué précédemment ) à une échelle$a=64$, et $\sigma=5$:
$$ \psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2} $$
est intégré et normalisé L1 : - (voir avertissement2 ci-dessous)
$$ \psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3} $$
La norme est décalée par un facteur scalaire, mais ce qui compte c'est la même chose pour différents$a$; vous pouvez vérifier que c'est le cas (pour les grands$a$, augmenter les limites d'intégration). Notez que c'est l' ondelette normalisée L1, pas l'ondelette intégrée$\psi_\text{int}(t)$; la distinction est essentielle, comme nous le verrons.
$\psi$ est ensuite convoluée avec le signal $s$, et aux différences finies et niées; d' ici , c'est:
$$ -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4} $$
En inspectant de près, cela équivaut à annuler l'ancienne intégration, avec certaines différences numériques (abordées plus tard), nous nous retrouvons donc avec
$$ \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5} $$
Enfin * sqrt(scale):
$$ \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6} $$
qui est donc finalement normalisé L2 . Nous pouvons le vérifier sur un signal de test; L2 est polarisé dans la représentation d'amplitude / puissance vers les basses fréquences, supprimant les hautes fréquences. Si ci-dessus est correct, non seulement pouvons-nous voir cela, mais appliquer un correctif en supprimant * sqrt(scale)du code:
Les mises en garde spécifiques à la mise en œuvre sont discutées ensuite.
Quelle est la fiabilité de la normalisation à travers les échelles? - Pas génial (- code2 ):
À droite, nous voyons le meilleur cas utilisant un nombre limité d'échantillons, à gauche, nous voyons l'implémentation du rééchantillonnage. Clairement en raison de la résolution intrinsèquement médiocre d'un petit nombre d'échantillons, nos normes ne fonctionnent pas bien pour les petites échelles - mais la variante recalculée maintient au moins une norme moyenne à peu près constante entre elles.
Plus intéressant encore, celui rééchantillonné apparaît soudainement quadratique après $\log_2(a)=6$; quel est le problème? Rappelez-vous "Cela commence ... à grande échelle$a=64$"; c'est quand nous rééchantillonnons avec plus d'échantillons que dans l'original ( int_psi):
à qui diffest
- Dans 1LR, nous voyons
int_psiàscale=256rééchantillonner les mêmes valeurs àint_psiplusieurs reprises, ce qui se manifeste par des zéros dansdiff(2LR). - Ensuite, nous mettons les valeurs au carré pour avoir une idée de la norme L2, et comparons avec les valeurs normalisées L2 dans 3LR; nous trouvons des valeurs normées soufflant des valeurs non formées hors de l'eau. Est-ce excessif ou approprié?
- Comparer avec recalculé à la même échelle; dans 4LR, on voit au moins la forme recalculée - ce n'est plus négligeable. Cela suggère que le rééchantillonnage surestime l'amplitude des ondelettes, et comme dans le calcul de la norme, cela est significativement décalé de la moitié des valeurs étant nulles en raison de
diff, alors la surestimation doit être grande.
Comme c'est génial? Nous pouvons comparer les amplitudes d'ondelettes de crête à travers les échelles des rééchantillonnages diff(int_psi_256)et recalculés non formés w( morlet2/scale); normer les changements à la fois par le même facteur ( *sqrt(scale)), ne fera aucune différence.
Assez sûrement, encore une fois à $a=64$nous voyons un comportement critique; l'amplitude de crête des saturés rééchantillonnés au minimum, tandis que la valeur recalculée continue de diminuer. La différence à$a=256$est plus que suffisant (~ 12db = 4x, si carré -> 16x) pour décaler les zéros diff, d'où la croissance de la norme L2.
Mais notez que la norme L1 reste bien comportée à des échelles plus élevées. Pour savoir pourquoi, il faut comparer précisément la proportion de zéros avec le facteur de surestimation. Nous savons déjà que ce dernier est 4x; pour les premiers, j'ai pris une région non plate (sinon c'est des zéros quand même) et j'ai trouvé le rapport des zéros au nombre d'échantillons: 0,75 . Ce qui signifie que 1 valeur sur 4 est non nulle, et chacune est 4x plus grande ... la proportion exacte nécessaire pour conserver la norme L1.
Cette proportion de zéros n'est pas très claire à partir des graphiques ci-dessus - donc à la place, les points:
Si diffdéfait cumsum, pourquoi faire l'un ou l'autre?
Bonne question. En fait, en me débarrassant des deux, j'ai trouvé que l'effet de dégradation de la qualité (fermeture à glissière) disparaît:
Mais quelque chose d'autre disparaît également - les échelles de haute fréquence (à gauche) et le signe se retourne (j'ai annulé cela en supprimant moins à coef = -). Je compense en partie cela en remplaçant * sqrt(scale)par 1 / sqrt(scale)(à droite), mais cela nous amène au cas normé L2 ( 1 / scalec'est exagéré), et c'est une modification sans fondement à la lumière de tout ce qui a été trouvé jusqu'à présent. Alors, est-ce possible? Bien sûr, on peut mettre en œuvre sans la int-diffdanse, mais je ne vais pas l'explorer ici.
Attention: norme effective
Rappel, convest ajusté à len(x), donc l' ondelette entière peut ne pas participer au calcul, alors que les calculs ci-dessus supposent le contraire. L'hypothèse est d'autant plus vraie que la longueur xest longue - mais voici ce qui se passe si len(x) == 512, avec la même ondelette precision=10(2 ^ 10 = 1024 int_psi), comparée à une ondelette recalculée:
Un gâchis! Les deux se déforment considérablement, ce qui implique peut-être qu'il est inhérent à des échelles suffisamment grandes épinglées contre des signaux courts. Une solution de contournement simple est de ne pas remplir le signal, mais je ne me porterais pas garant sans comprendre les mises en garde (ce que je ne fais pas pour les ondelettes, mais pour DFT ).
De plus, notez que la norme n'est pas seulement par échelle, mais aussi par équipe ($b$); moins de l'ondelette est dans le signal quand il vient juste d'entrer par opposition à quand son centre est au centre du signal.
Mais en le gardant par échelle, j'avoue avoir triché dans l'exemple précédent avec $f=4$; au départ j'ai choisi$f=1$, mais l'intrigue de la norme L1 n'avait pas l'air sympa, et je ne comprenais pas pourquoi. Depuis que je le fais, voici à la fois L1 et L2 par échelle pour$f=1,64$:
Comparez avec l'intrigue précédente; notez comment les modes de gauche (échelles basses) des deux normes ne sont pas affectés, seuls les modes haute (basse fréquence) sont maintenant beaucoup plus grands (L2) et plus bas (L1), comme prévu.
Attention 2: * stepne fait pas la normalisation L1
Au lieu de cela, il résulte de diffl'ingestion de l'ondelette intégrée rééchantillonnée; le * stepseul définit l'échelle initiale correcte et peut être omis sans perdre la norme. Cela ne se produit pas avec le recalcul de l'ondelette ou le rééchantillonnage de l'ondelette non intégrée.
(Notez que recalculé scipy.signal.morlet2est normalisé en L2, mais uniquement parce qu'il est explicitement normalisé; le diff-réchantillon- intest "auto-normalisé").