PyWavelets CWT: нормализация? Против Scipy?
Связанные . Реализуемое уравнение нормализуется 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) $$
так почему pywt._cwtже * sqrt (scale) ? Кроме того, как это нормализовать (L1, L2 и т. Д.)?
Изменить : просмотрев cwt Scipy , я понял, что он реализует именно то, что я назвал «идеальным» в ответе, таким образом, это то же самое, что и сравнение двух (за исключением нескольких очевидных деталей). Я не знаю, в каких сценариях алгоритм pywt лучше, но в целом - легкая победа для Scipy . Я открыл PR, чтобы устранить ключевой недостаток pywt.
Edit2 : обнаружил явное преимущество pywt над scipy в наименее ожидаемом месте: низкие масштабы . Адресовано здесь .
Редактирование 3 : после закрепления минимума precisionв pywt, я фактически нашел pywtвыше scipyво всех из множества синтетических сигналов, но squeezepy«ы всегда бить обоих. Тесты здесь , скоро будут.
Ответы
Давайте проследим математику от инкубации до доставки. Он начинается с psiмасштабирования morlet2(как показано ранее ) в масштабе$a=64$, а также $\sigma=5$:
$$ \psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2} $$
получает интегрированный и L1-нормированы : - (см caveat2 ниже)
$$ \psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3} $$
Норма выключена с помощью скалярного множителя, но важно то , что это же для разных$a$; вы можете проверить это так (для больших$a$, увеличьте пределы интеграции). Обратите внимание, что это вейвлет , нормализованный по L1, а не интегрированный вейвлет.$\psi_\text{int}(t)$; различие существенно, как мы увидим.
$\psi$ затем сворачивается с сигналом $s$, и конечно-разностные и отрицательные; от здесь , это:
$$ -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4} $$
Если присмотреться, это эквивалентно отмене прежней интеграции с некоторыми числовыми различиями (рассмотренными позже), поэтому мы остаемся с
$$ \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5} $$
И наконец * sqrt(scale):
$$ \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6} $$
который, таким образом, в конечном итоге L2-нормализован . Мы можем проверить это на тестовом сигнале; L2 смещен в представлении амплитуды / мощности в сторону более низких частот, подавляя высокие частоты. Если это верно, мы не только можем это увидеть, но и применить исправление, удалив * sqrt(scale)из кода:
Предостережения, связанные с реализацией, обсуждаются далее.
Насколько надежна нормализация по шкалам? - Не отлично (- code2 ):
Справа мы видим лучший случай с использованием ограниченного количества выборок, слева мы видим реализацию передискретизации. Ясно, что из-за изначально плохого разрешения небольшого количества выборок наши нормы не подходят для низких шкал, но пересчитанный вариант по крайней мере поддерживает среднюю норму примерно постоянной для всех них.
Что еще более интересно, передискретизированное изображение внезапно становится квадратичным после $\log_2(a)=6$; в чем дело? Напомним: «Все начинается ... в масштабе$a=64$"; это когда мы повторно сэмплируем больше сэмплов, чем в original ( int_psi):
которого diffявляется
- В 1LR мы видим, что
int_psiatscale=256- это повторная выборка одних и тех же значенийint_psiнесколько раз, что проявляется в виде нулей вdiff(2LR). - Затем мы возводим значения в квадрат, чтобы получить представление о норме L2, и сравниваем с L2-нормализованными значениями в 3LR; мы обнаруживаем, что нормированные значения вытесняют ненормированные значения из воды. Это чрезмерно или уместно?
- Сравнить с пересчитанными в том же масштабе; в 4LR мы, по крайней мере, видим пересчитанную форму - ею больше нельзя пренебрегать. Это предполагает, что передискретизация переоценивает амплитуду вейвлета, и поскольку при вычислении нормы это значительно компенсируется половиной значений, равных нулю, из-за того
diff, что переоценка должна быть большой.
Как здорово? Мы можем сравнивать пиковые амплитуды вейвлетов по масштабам без нормированной передискретизации diff(int_psi_256)и пересчета w( morlet2/scale); нормализация изменений обоих по одному и тому же фактору ( *sqrt(scale)) не будет иметь значения.
Конечно же, снова в $a=64$мы видим критическое поведение; пиковая амплитуда передискретизированного насыщается до минимума, тогда как пересчитанная продолжает уменьшаться. Разница на$a=256$более чем достаточно (~ 12 дБ = 4x, если возведено в квадрат -> 16x), чтобы компенсировать нули diff, поэтому норма L2 растет.
Но обратите внимание, что норма L1 остается хорошей на более высоких масштабах. Чтобы выяснить, почему, мы должны точно сравнить долю нулей с коэффициентом завышения. Мы уже знаем, что последнее - 4х; для первого я взял неплоскую область (иначе все равно нули) и нашел отношение нулей к количеству отсчетов: 0,75 . Это означает, что 1 из 4 значений не равны нулю, и каждое в 4 раза больше ... точная пропорция, необходимая для сохранения нормы L1.
Эта пропорция нулей не очень ясна из приведенных выше графиков - поэтому вместо этого точки:
Если diffотменяет cumsum, зачем вообще то и другое?
Хороший вопрос. Фактически, избавившись от обоих, я обнаружил, что эффект ухудшения качества (застежка-молния) исчез:
Но пропадает еще кое-что - шкала высоких частот (слева), и знак переворачивается (я это убрал, убрав минус у coef = -). Я частично компенсирую это, заменяя его * sqrt(scale)на 1 / sqrt(scale)(справа), но это подводит нас к нормированному случаю L2 ( 1 / scaleэто избыточно) и является необоснованной модификацией в свете всего, что было найдено до сих пор. Так можно ли это сделать? Конечно, можно реализовать без int-diffтанца, но я не буду здесь его исследовать.
Предупреждение: эффективная норма
Напомним, convобрезается до len(x), поэтому весь вейвлет может не участвовать в вычислениях, тогда как вышеприведенные вычисления предполагают иное. Предположение тем вернее, чем длиннее x- но вот что произойдет, если len(x) == 512с тем же вейвлетом precision=10(2 ^ 10 = 1024 int_psi) по сравнению с пересчитанным вейвлетом:
Беспорядок! Оба искажения значительно искажаются, возможно, подразумевая, что это присуще достаточно большим масштабам, привязанным к коротким сигналам. Простой обходной путь - заполнение сигнала нулями, но я бы не поручился за него, не понимая предостережений (которые я не делаю для вейвлетов, но делаю для DFT ).
Кроме того, обратите внимание, что норма не только по шкале, но и по смене ($b$); меньше вейвлета находится в сигнале, когда он только входит, в отличие от того, когда его центр находится в центре сигнала.
Но, сохраняя масштаб, признаюсь, в предыдущем примере я обманул $f=4$; изначально я выбрал$f=1$, но график нормы L1 выглядел некрасиво, и я не понимал почему. Так как сейчас я делаю, вот L1 и L2 для каждой шкалы для$f=1,64$:
Сравните с более ранним сюжетом; обратите внимание, что левые моды (низкие шкалы) обеих норм не затронуты, только высокие шкалы (низкие частоты) теперь намного больше (L2) и ниже (L1), как и ожидалось.
Предостережение 2: * stepне выполняет L1-нормализацию
Вместо этого он является результатом diffповторной дискретизации интегрированного вейвлета; * stepтолько устанавливает правильную начальную шкалу, и может быть опущена без потери нормы. Этого не происходит при повторном вычислении вейвлета или повторной дискретизации неинтегрированного вейвлета.
(Обратите внимание, что перерасчет scipy.signal.morlet2является L2-нормализованным, но только потому, что он явно нормализован; параметр diff-resample- int«самонормализован»).