Реализация PyWavelets CWT

Question

Я пытаюсь понять реализацию PyWavelets непрерывного вейвлет-преобразования и ее сравнение с более «базовой» версией, которую я закодировал и предоставил здесь. Особенно:

1. Как интегрированный вейвлет, а затем разностный пост-свертка эквивалентен свертке без int + diff?
2. Как передискретизация интегрированного вейвлета при увеличении разрешения эквивалентна изменению его масштаба?

---

Наивная реализация :

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$синусоиды :

Answer 1

Сначала см. «Наивный анализ» в разделе ниже. На PyWavelets: на Github был обнаружен алгоритм, основанный на старой реализации MATLAB , но он не предоставляет подробностей о кодировании вейвлет-ядра. Итак, я принялся исследовать; различные визуальные эффекты и код аннотаций здесь . Ниже приводятся основные выводы. - Весь код

См. «Разбор PyWavelets» ниже. Он реализует следующее, а именно уравнение 4:

---

(1) Как интегрированный вейвлет, а затем разностная свертка после свертки эквивалентны свертке без int + diff?

Формулировка находится справа, получая из левого уравнения; Вместо того, чтобы постоянно интегрировать все$\mathcal{R}$, интеграл разбивается на отрезки s, суммируя $k$. Вход считается постоянным в течение$[k, k+1]$ (хотя не должно быть $)$?), поскольку он дискретный, тогда как вейвлет остается непрерывным. Это позволяет вывести вход за пределы интеграла; визуально:

Итак, как нам это кодировать? Из уравнения 4 ясно, что diffгде- то нам нужен . pywtделает np.diff(conv), что в основном np.diff(coef); почему в Наиве вроде бы берется разница в коэффициенте мнимого?

В уравнении 4 у нас есть интегрированный вейвлет до k + 1минус до k, умноженный на вейвлет. Это эквивалентно произведению с вейвлетом на b - 1и b, соответственно (подключите t=k+1и t=k, сравните с подключением b=b-1и b=bвместо`).

Ответ тогда в conv; предположим, что он уже обрезан, так что это только тогда, когда вейвлет находится по крайней мере на полпути внутри сигнала. Тогда каждая точка в convявляется произведением вейвлета с частью входа, в другом tau(точнее b). Итак, две соседние точки conv- это произведения с вейвлетами на bи 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} $$

Тем не менее, мы ищем минус между t=k+1и tили большее tи меньшее t, что является минусом вышеизложенного, что объясняет шаг 7 в PW Breakdown.

---

(2) Как передискретизация интегрированного вейвлета с увеличением разрешения эквивалентна изменению его масштаба?

См. Рисунок под (5) в PW Breakdown. pywtпохоже, "масштабируется" через количество выборок, которые определяют вейвлет, а не пересчитывает вейвлет с расширенным аргументом. ... хотя, визуально, это передискретизация явно делает растянуть вейвлет, растягивая его по отношению к входу.

Так что это действительно форма масштабирования; но как это сравнить с наивным? Naive пересчитывается более интуитивно по уравнениям, используя фиксированную длину вейвлета и передавая масштабированный аргумент. Тогда что лучше? На противоположный вопрос ответить яснее; подумайте, что происходит на более высоких и низших уровнях:

1. Наивно выше : вейвлет слишком сильно расширяется, покидая рамку выборки намного раньше, чем ее нулевые хвосты.
2. pywt выше : вейвлет вполне подойдет ; проблема с conv. Длина вейвлета зафиксирована на уровне 1024, поэтому, если входной сигнал короче, вейвлеты более высокого масштаба никогда не смогут полностью умножить сигнал. Чем больше несоответствие, тем больше вейвлет "видится" сигналом, как "Наивно выше"; это видно на тепловых картах вопроса, отличающихся вертикальным смещением.
   • Максимальный масштаб вейвлета также не ограничен; len(j)растет пропорционально масштабу (но есть проблема нормализации, см. здесь ).
3. Наивный нижний : все хорошо.
4. pywt lower : вейвлет имеет маленькое разрешение.
5. Фиксированная и переменная длина окна : более высокие масштабы должны соответствовать большему расширению или отличному от нуля вейвлету относительно входа . Оба достигают этого, за исключением того, что Naive резко обнуляется на более высоких масштабах.

В обоих случаях проблемы более высокого масштаба решаются путем ограничения максимального масштаба; если оба находятся в своих полностью безопасных зонах, наивная реализация кажется предпочтительной из-за более высокого разрешения. Однако эта «безопасная зона» для Naive непомерно мала, если считать «безопасным» отсутствие прыжка к нулю на хвосте.

Напротив, pywtпроблему легче решить, увеличив максимальную длину вейвлета, что было предложено для избавления от артефактов, похожих на застежки-молнии (которые можно увидеть на графиках вопроса). Это, однако, ухудшило бы производительность в более высоком масштабе для коротких входов - но на практике оно вряд ли будет меньше 1024.

---

PyWavelets разбивка :

1. Вейвлет до интеграции точно совпадает с показанным фрагментом кода, который является приближением полного реального Морлета (используемого Naive), предполагая $\sigma > 5$в Вики .
2. pywtинтегрирует реальный переходник Морле np.cumsum(psi) * step, учитывая разный размер шага
3. Интегрированный вейвлет int_psiповторно используется для всехscales
4. Для каждого scale, так же int_psiкак понижается при увеличении разрешения, причем , jкак целое число линейно между различной же мин и макс (почти) для всех scales:

1. Результирующий вейвлет для увеличения scale, бок о бок с наивным (на будущее):

1. Результат свертки сначала diffd,
2. отрицается,
3. умноженный на sqrt(scale),
4. затем обрезали с обеих сторон , чтобы либо (а) включают в себя только части , где вейвлет по меньшей мере , на полпути «внутри» сигнала, или (б) заставить его длины к len(x).

---

Наивная поломка : если вы не знакомы с CWT / STFT, я настоятельно рекомендую здесь части 1-3 .

Нам понадобится (1) тип вейвлета; (2) длина вейвлета; (3) приращение / перекрытие вейвлетов. Мы будем использовать настоящий Морле , 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(«приращение») во внешнем цикле, там масштабирование, но они эквивалентны.

Наконец, формы при тусклом сдвиге во времени различаются; реализация вопроса обрезается так, что включаются только те части свертки, где вейвлет находится по крайней мере на полпути "внутри входа"; этот урезает дальше до того места, где включен полный вейвлет.

---

НЕРЕШЕННЫЕ :

1. Почему * sqrt(scale)(шаг 8 pywt)? Уравнение 4 показывает 1 / sqrt(scale). Я рассмотрел * scaleправило для каждой цепочки, за исключением того, что мы интегрировали вейвлет численно ... Имеет смысл иметь индивидуальную обработку; мы получаем далеко с не интегрировании за масштаба вейвлета , потому что масштаб за вейвлет ... тот же передискретизируется Вейвлет? Решено .