การใช้งาน 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" ด้านล่าง กำลังดำเนินการต่อไปนี้ - คือ Eq 4:

---

(1) เวฟเล็ตแบบบูรณาการเป็นอย่างไรแล้วโพสต์คอนเวอลั่มที่แตกต่างกันอย่างไรเทียบเท่ากับการชักโดยไม่มี int + diff

สูตรอยู่ทางขวามาจากสมการทางซ้าย แทนที่จะรวมอย่างต่อเนื่องทั้งหมด$\mathcal{R}$อินทิกรัลถูกแบ่งออกเป็นเซกเมนต์ ssumming มากกว่า $k$. อินพุตจะถือว่าคงที่มากกว่า$[k, k+1]$ (แม้ว่าจะไม่ควรก็ตาม $)$?) เนื่องจากไม่ต่อเนื่องในขณะที่เวฟเล็ตจะยังคงต่อเนื่อง สิ่งนี้ทำให้สามารถย้ายอินพุตออกนอกอินทิกรัล สายตา:

แล้วเราจะเขียนโค้ดนี้อย่างไร? จาก Eq 4 ชัดเจนว่าเราต้องการdiffที่ไหนสักแห่ง pywtไม่np.diff(conv)ซึ่งโดยพื้นฐานnp.diff(coef)แล้ว ทำไมค่าสัมประสิทธิ์ที่จะเป็นจริงใน Naive จึงแตกต่างกัน

ใน Eq 4 เรามีเวฟเล็ตแบบบูรณาการ up-to k + 1, ลบ up-to 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 ต่ำกว่า : เวฟเล็ตมีความละเอียดน้อย
5. คงเทียบกับระยะเวลาในหน้าต่างตัวแปร : เครื่องชั่งน้ำหนักที่สูงขึ้นจะต้องสอดคล้องกับการขยายตัวมากขึ้นหรือเป็น-ภัณฑ์ของเวฟเทียบกับการป้อนข้อมูล ทั้งสองทำสิ่งนี้สำเร็จยกเว้น Naive หยุดเป็นศูนย์ทันทีที่สเกลที่สูงขึ้น

ปัญหาระดับที่สูงขึ้นในทั้งสองได้รับการแก้ไขโดยการ จำกัด มาตราส่วนสูงสุด หากทั้งสองอยู่ในเขตปลอดภัยอย่างสมบูรณ์การใช้งานแบบไร้เดียงสาดูเหมือนจะดีกว่าเนื่องจากมีความละเอียดสูงกว่า อย่างไรก็ตาม "เขตปลอดภัย" สำหรับ Naive นี้มีขนาดเล็กมากหากเราพิจารณาว่า "ปลอดภัย" ไม่ใช่กระโดดไปที่ศูนย์ที่หาง

ในทางตรงกันข้ามpywtปัญหาจะแก้ไขได้ง่ายขึ้นโดยการเพิ่มความยาวเวฟเล็ตสูงสุดซึ่งแนะนำให้กำจัดสิ่งประดิษฐ์ที่มีลักษณะคล้ายซิป (ซึ่งสามารถดูได้ในพล็อตคำถาม) อย่างไรก็ตามนั่นทำให้ประสิทธิภาพที่สูงขึ้นสำหรับอินพุตสั้น ๆ แย่ลง แต่ในทางปฏิบัติไม่น่าจะสั้นกว่า 1024

---

รายละเอียด PyWavelets :

1. Wavelet ก่อนการรวมจะตรงกับโค้ดหยดที่แสดงซึ่งเป็นการประมาณค่า Morlet จริงที่สมบูรณ์ (ใช้โดย Naive) โดยสมมติว่า $\sigma > 5$ในวิกิพีเดีย
2. pywtรวม Morlet จริงผ่านการnp.cumsum(psi) * stepบัญชีสำหรับขนาดขั้นตอนที่แตกต่างกัน
3. เวฟเล็ตint_psiในตัวถูกนำกลับมาใช้ใหม่สำหรับทุกคนscales
4. สำหรับแต่ละscaleรายการint_psiจะมีการสุ่มตัวอย่างแบบเดียวกันที่ความละเอียดที่เพิ่มขึ้นโดยjจำนวนเต็มจะแปรผันตามเส้นตรงระหว่างค่าต่ำสุดและสูงสุดเท่ากัน (เกือบ) สำหรับทั้งหมดscales:

1. เวฟเล็ตที่เป็นผลลัพธ์สำหรับการเพิ่มขึ้นscaleเคียงข้างกันอย่างไร้เดียงสา (ในภายหลัง):

1. ผลการแปลงครั้งแรกdiff'd,
2. เมื่อตะกี้
3. คูณด้วยsqrt(scale),
4. แล้วตัดแต่งทั้งสองด้านอย่างใดอย่างหนึ่ง (ก) เท่านั้นรวมถึงชิ้นส่วนที่เวฟเป็นอย่างน้อยครึ่งหนึ่ง "ภายในสัญญาณ" หรือ (ข) len(x)บังคับให้ระยะเวลาในการ

---

พังทลายหน่อมแน้ม : ถ้าไม่คุ้นเคยกับ CWT / STFT ผมขอแนะนำส่วนที่ 1-3 ที่นี่

เราต้องการ (1) ประเภทเวฟ; (2) ความยาวเวฟ; (3) การเพิ่มขึ้นของเวฟเล็ต / การทับซ้อนกัน เราจะใช้Morletจริง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("Increment") ในวงนอกมีมาตราส่วน แต่ทั้งสองเท่ากัน

สุดท้ายรูปร่างตามไทม์ชิฟสลัวแตกต่างกัน การดำเนินการของคำถามจดจ้องเพื่อให้มีการรวมเฉพาะบางส่วนของคอนโวลูชั่นที่เวฟเล็ตอย่างน้อยครึ่งหนึ่ง "อยู่ในอินพุต" อันนี้จดจ้องไปที่ที่รวมเวฟเล็ตแบบเต็ม

---

ยังไม่ได้แก้ไข :

1. ทำไม* sqrt(scale)(ขั้นตอนที่ 8 pywt)? สมการที่ 4 1 / sqrt(scale)แสดงให้เห็นว่า ฉันพิจารณา* scaleตามกฎลูกโซ่ยกเว้นว่าเราได้รวมเวฟเล็ตเป็นตัวเลข ... สมเหตุสมผลที่จะมีการรักษาต่อสเกล เราจะหลีกเลี่ยงการไม่รวมเวฟเล็ตต่อสเกลได้หรือไม่เพราะเวฟเล็ตต่อสเกลคือ ... เวฟเล็ตที่สุ่มตัวอย่างเหมือนกัน? แก้ไขแล้ว