Triển khai PyWavelets CWT
Tôi tìm cách hiểu việc triển khai Biến đổi Wavelet Liên tục của PyWavelets và cách nó so sánh với phiên bản 'cơ bản' hơn mà tôi đã mã hóa và cung cấp ở đây. Đặc biệt:
- Làm thế nào được tích hợp wavelet, sau đó biến đổi sau chuyển đổi, tương đương với chuyển đổi không có int + diff?
- Làm thế nào để lấy lại mẫu wavelet tích hợp ở độ phân giải tăng lên tương đương với việc thay đổi tỷ lệ của nó?
Thực hiện ngây thơ :
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 (được cắt vào các phần có liên quan) :
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
So sánh trên $f=1, 4$hình sin :
Trả lời
Trước tiên, hãy xem "Phá vỡ Naive" trong phần bên dưới. Đối với PyWavelets: thuật toán được tìm thấy trên Github bắt nguồn từ việc triển khai MATLAB cũ , nhưng nó không cung cấp chi tiết về việc mã hóa hạt nhân wavelet. Vì vậy, tôi bắt đầu khám phá; nhiều hình ảnh và chú thích mã ở đây . Những phát hiện chính sau đây. - Tất cả mã
Xem "Phân tích PyWavelets" bên dưới. Nó đang thực hiện những điều sau - cụ thể là Công thức 4:
(1) Làm thế nào là Wavelet tích hợp, sau đó phân biệt sau chuyển đổi, tương đương với chuyển đổi không có int + diff?
Công thức ở bên phải, suy ra từ bên trái, phương trình; Thay vì tích hợp liên tục trên tất cả$\mathcal{R}$, tích phân được chia thành các phân đoạn gộp lại $k$. Đầu vào được giả định là không đổi$[k, k+1]$ (mặc dù không nên $)$?), vì nó rời rạc, trong khi wavelet được giữ liên tục. Điều này cho phép di chuyển đầu vào bên ngoài tích phân; trực quan:
Vì vậy, làm thế nào để chúng tôi mã này? Rõ ràng từ phương trình 4 rằng chúng ta cần một diffnơi nào đó. pywthiện np.diff(conv), về cơ bản là np.diff(coef); tại sao nó dường như có sự khác biệt của hệ số sẽ có trong Naive?
Trong phương trình 4, chúng ta có Wavelet tích hợp tăng lên k + 1, trừ lên đến k, nhân với wavelet. Điều này tương đương với sản phẩm có wavelet tại b - 1và b, tương ứng (cắm vào t=k+1và t=k, so sánh với cắm vào b=b-1và b=bthay vào đó ').
Câu trả lời sau đó là trong conv; giả sử nó đã được cắt bớt để chỉ khi wavelet nằm trong ít nhất một nửa bên trong tín hiệu. Sau đó, mọi điểm trong convlà sản phẩm của wavelet với một phần đầu vào, ở một điểm khác tau(hay đúng hơn b). Vì vậy, hai điểm liền kề trong convlà các sản phẩm có wavelets tại bvà 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} $$
Tuy nhiên, chúng tôi tìm kiếm trừ giữa t=k+1và t, hoặc lớn hơn tvà thấp hơn t, là phủ định của ở trên, do đó giải thích bước 7 trong Phân tích PW.
(2) Việc lấy lại mẫu wavelet tích hợp ở độ phân giải tăng tương đương với việc thay đổi tỷ lệ của nó như thế nào?
Xem hình dưới (5) trong Phân tích PW. pywtdường như "chia tỷ lệ" thông qua số lượng mẫu xác định wavelet, thay vì tính toán lại wavelet với một đối số giãn ra. ... tuy nhiên, về mặt trực quan, việc lấy lại mẫu này rõ ràng đã làm giãn wavelet, kéo dài nó so với đầu vào.
Vì vậy, nó thực sự là một hình thức mở rộng quy mô; nhưng làm thế nào để nó so sánh với ngây thơ? Naive được tính toán lại trực quan hơn cho mỗi phương trình, sử dụng độ dài wavelet cố định và truyền đối số được chia tỷ lệ. Sau đó, cái nào tốt hơn? Câu hỏi ngược lại rõ ràng hơn để trả lời; xem xét những gì xảy ra ở quy mô cao hơn và thấp hơn:
- Ngây thơ cao hơn : wavelet giãn ra quá nhiều, thoát ra khỏi khung lấy mẫu trước khi có đuôi 0 của nó.
- pywt cao hơn : wavelet là tốt; vấn đề là với
conv. Độ dài Wavelet được cố định ở 1024, vì vậy nếu đầu vào ngắn hơn bất kỳ, thì các wavelet có quy mô cao hơn không bao giờ có thể nhân đầy đủ tín hiệu. Chênh lệch càng lớn, wavelet càng được "nhìn thấy" tương tự như "Naive cao hơn" bởi tín hiệu; điều này có thể được nhìn thấy trong bản đồ nhiệt của câu hỏi khác nhau theo sự dịch chuyển dọc.- Quy mô tối đa của wavelet cũng không giới hạn;
len(j)tăng tỷ lệ thuận với quy mô (nhưng có một vấn đề về chuẩn hóa, xem tại đây ).
- Quy mô tối đa của wavelet cũng không giới hạn;
- Ngây thơ hạ : tất cả đều tốt.
- pywt thấp hơn : wavelet có độ phân giải nhỏ.
- Độ dài cửa sổ cố định và thay đổi : tỷ lệ cao hơn phải tương ứng với độ giãn nở lớn hơn, hoặc không phải là khác, của wavelet so với đầu vào . Cả hai đều thực hiện được điều này, ngoại trừ Naive đột ngột dừng về 0 ở các quy mô cao hơn.
Các vấn đề về quy mô cao hơn trong cả hai đều được giải quyết bằng cách giới hạn quy mô tối đa; nếu cả hai đều ở trong vùng hoàn toàn an toàn của chúng, việc triển khai Naive có vẻ thích hợp hơn do độ phân giải cao hơn. Tuy nhiên, "vùng an toàn" này đối với Naive là rất nhỏ, nếu chúng ta coi "an toàn" là không nhảy về số 0 ở đuôi.
Ngược lại, pywtvấn đề của nó được giải quyết dễ dàng hơn bằng cách tăng chiều dài wavelet tối đa, điều này được đề xuất để loại bỏ các hiện vật giống dây kéo (có thể thấy trong các ô của câu hỏi). Tuy nhiên, điều đó sẽ làm xấu đi hiệu suất ở quy mô cao hơn đối với các đầu vào ngắn - nhưng trên thực tế, nó không chắc sẽ ngắn hơn 1024.
Phân tích PyWavelets :
- Wavelet, trước khi tích hợp, khớp chính xác với khối mã được hiển thị, là một xấp xỉ của Morlet thực hoàn chỉnh (được sử dụng bởi Naive) $\sigma > 5$trong Wiki .
pywttích hợp Morlet thực thông quanp.cumsum(psi) * step, tính toán kích thước bước khác biệt- Wavelet tích hợp
int_psi, được sử dụng lại cho tất cảscales - Đối với mỗi cái
scale, cái giống nhauint_psiđược lấy mẫu lại ở độ phân giải tăng dần, vớijdạng số nguyên thay đổi tuyến tính giữa cùng một mức tối thiểu và tối đa (gần như) cho tất cảscales:
- Wavelet kết quả để tăng
scale, song song với ngây thơ (cho phần sau):
- Kết quả tích chập đầu tiên
diff'd, - phủ định,
- nhân với
sqrt(scale), - sau đó được cắt ở cả hai bên để (a) chỉ bao gồm các phần mà wavelet nằm ít nhất một nửa "bên trong tín hiệu", hoặc (b) buộc chiều dài của nó lên
len(x).
Naive Breakdown : nếu không quen với CWT / STFT, tôi thực sự giới thiệu phần 1-3 ở đây .
Chúng ta cần (1) loại wavelet; (2) chiều dài wavelet; (3) gia tăng / chồng chéo wavelet. Chúng tôi sẽ sử dụng Morlet thực , 96 mẫu và tăng thêm 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
Hình ảnh hóa cho scale=2:
Sự khác biệt giữa cách triển khai này và cách triển khai ngây thơ trong câu hỏi là, cách sử dụng sau này win_inc=1, tương đương với việc xoay vòng với một wavelet lật và cắt bớt kết quả để chỉ bao gồm các phần mà wavelet hoàn toàn "bên trong tín hiệu". Ngoài ra, ở đây chúng ta có tau("gia số") trong vòng lặp bên ngoài, có quy mô, nhưng hai tương đương.
Cuối cùng, hình dạng dọc theo thời gian mờ khác nhau; phần thực hiện của câu hỏi cắt bớt để chỉ các phần của tích chập trong đó wavelet nằm ít nhất nửa "bên trong đầu vào" mới được đưa vào; cái này cắt tỉa thêm đến nơi có đầy đủ wavelet.
CHƯA GIẢI QUYẾT :
Tại saoĐã giải quyết .* sqrt(scale)(bước 8 pywt)? Phương trình 4 cho thấy1 / sqrt(scale). Tôi đã xem xét* scaletheo quy tắc chuỗi, ngoại trừ chúng tôi đã tích hợp wavelet bằng số ... Có ý nghĩa để xử lý theo quy mô; chúng ta có tránh được việc không tích hợp wavelet theo tỉ lệ không vì wavelet theo tỉ lệ là ... cùng một wavelet được lấy mẫu lại?