📋 Week 3: Properties of the Continuous Wavelet Transform 📝 Outline Goals

Inversion (Reconstruction) Formula

Energy Conservation (Parseval’s Theorem for CWT)

Scalogram Interpretation

Worked Example

Exercises

Energy Estimation

Energy of Continuous Wavelet

Scalogram Vs Spectrogram

Key takeaway

• Scalogram ≠ Spectrogram, though both show when and at what “frequency-like” band energy occurs.

• Scalograms use wavelets to give a zoomed-in view at high frequencies (good time detail) and a zoomed-out view at low frequencies (good frequency detail).

• Spectrograms slice the signal into equal-length windows, giving constant resolution everywhere.

---

Here’s a self‐contained script that you can drop into your venv (with PyWavelets ≥ 1.1.0) and run end-to-end. It will:

# continuous_wavelet_analysis.py

import numpy as np
import matplotlib.pyplot as plt
import pywt
from scipy.signal import chirp
from scipy.interpolate import interp1d

# 1. Generate a non-stationary signal: two chirps + noise
fs = 1000            # Hz
T  = 1.0             # seconds
t  = np.linspace(0, T, int(fs*T), endpoint=False)

signal = np.zeros_like(t)
#  up-chirp from 50→200 Hz
signal[:fs//2] = chirp(t[:fs//2], f0=50,  f1=200, t1=T/2, method='linear')
# down-chirp from 300→100 Hz
signal[fs//2:] = chirp(t[fs//2:], f0=300, f1=100, t1=T/2, method='linear')
# add white noise
signal += 0.2 * np.random.randn(len(t))

# 2. Time-domain energy
E_time = np.trapz(np.abs(signal)**2, t)
print(f"Time-domain energy: {E_time:.6f}")

# 3. Continuous wavelet transform
wavelet = 'morl'
scales  = np.arange(1, 128)
coefs, freqs = pywt.cwt(signal, scales, wavelet,
                        sampling_period=1/fs)

# 4. Raw time–scale energy (no 1/Cψ normalization)
db = t[1] - t[0]
da = scales[1] - scales[0]
E_raw = np.sum(
    np.abs(coefs)**2 * (db * da / (scales[:, None]**2))
)
print(f"Raw time–scale energy (∫∫|W|²/(a²)db da): {E_raw:.6f}")

# 5. Estimate admissibility constant Cψ ≃ E_raw/E_time
Cpsi_est = E_raw / E_time
print(f"Estimated admissibility constant C_ψ ≃ {Cpsi_est:.6f}")

# 6. Verify normalized energy matches time-domain
E_norm = E_raw / Cpsi_est
print(f"Normalized time–scale energy: {E_norm:.6f}   <-- should ≃ Time-domain")

# 7. Inverse CWT (built-in if possible, else manual fallback)
def manual_icwt(coefs, scales, wname, Cpsi, t, db, da):
    wav = pywt.ContinuousWavelet(wname)
    # sample mother wavelet at scale=1
    psi_vals, psi_x = wav.wavefun(level=10)
    interp_psi = interp1d(psi_x, psi_vals,
                          bounds_error=False, fill_value=0.0)
    x_rec = np.zeros_like(t)
    for j, a in enumerate(scales):
        norm = a**1.5
        for i, b in enumerate(t):
            W = coefs[j, i]
            arg = (t - b) / a
            x_rec += W * (interp_psi(arg)/norm) * db * da
    return x_rec / Cpsi

if hasattr(pywt, 'icwt'):
    signal_rec = pywt.icwt(coefs, scales, wavelet,
                           sampling_period=1/fs)
else:
    print("⚠️ pywt.icwt not available, using slow manual inverse")
    signal_rec = manual_icwt(coefs, scales, wavelet,
                             Cpsi_est, t, db, da)

# 8. Reconstruction performance
mse = np.mean((signal - signal_rec)**2)
snr = 10 * np.log10(np.sum(signal**2) / np.sum((signal - signal_rec)**2))
print(f"Reconstruction MSE: {mse:.6e}")
print(f"Reconstruction SNR: {snr:.2f} dB")

# 9. Plot results
plt.figure(figsize=(10,3))
plt.plot(t, signal,     label='Original',    alpha=0.7)
plt.plot(t, signal_rec, label='Reconstructed',alpha=0.7)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original vs. Reconstructed Signal')
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()

Output