diff --git a/examples/common.py b/examples/common.py
new file mode 100644
index 0000000..30bfbd4
--- /dev/null
+++ b/examples/common.py
@@ -0,0 +1,98 @@
+from typing import Optional, Tuple, Union
+
+import bilby
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.interpolate import interp1d
+
+
+def compute_wavelet_snr(h: "Wavelet", PSD: "Wavelet") -> float:
+    """Compute the SNR of a model h[ti,fi] given data d[ti,fi] and PSD[ti,fi].
+
+    SNR(h) = Sum_{ti,fi} [ h_hat[ti,fi] d[ti,fi] / PSD[ti,fi]
+
+    Parameters
+    ----------
+    h : np.ndarray
+        The model in the wavelet domain (binned in [ti,fi]).
+    d : np.ndarray
+        The data in the wavelet domain (binned in [ti,fi]).
+    PSD : np.ndarray
+        The PSD in the wavelet domain (binned in [ti,fi]).
+
+    Returns
+    -------
+    float
+        The SNR of the model h given data d and PSD.
+
+    """
+    snr_sqrd = np.nansum((h * h) / PSD)
+    return np.sqrt(snr_sqrd)
+
+
+def compute_frequency_optimal_snr(h_freq, psd, duration) -> float:
+    """
+    A18 from Veitch et al. 2009
+    https://arxiv.org/abs/0911.3820
+    """
+    snr_sqrd = __noise_weighted_inner_product(
+        aa=h_freq, bb=h_freq, power_spectral_density=psd, duration=duration
+    ).real
+    return np.sqrt(snr_sqrd)
+
+
+def __noise_weighted_inner_product(aa, bb, power_spectral_density, duration):
+    integrand = np.conj(aa) * bb / power_spectral_density
+    return (4 / duration) * np.sum(integrand)
+
+
+def evolutionary_psd_from_stationary_psd(
+    psd: np.ndarray,
+    psd_f: np.ndarray,
+    f_grid,
+    t_grid,
+) -> "Wavelet":
+    """
+    PSD[ti,fi] = PSD[fi] * delta_f
+    """
+
+    Nt = len(t_grid)
+    delta_F = f_grid[1] - f_grid[0]
+    delta_T = t_grid[1] - t_grid[0]
+
+    freq_data = psd
+    nan_val = np.max(freq_data)
+    psd_grid = (
+        interp1d(
+            psd_f,
+            freq_data,
+            kind="nearest",
+            fill_value=nan_val,
+            bounds_error=False,
+        )(f_grid)
+        # * delta_F
+    )
+
+    # repeat the PSD for each time bin
+    psd_grid = np.repeat(psd_grid[None, :], Nt, axis=0)
+    return psd_grid
+
+
+def get_wavelet_bins(duration, data_len, Nf, Nt):
+    """Get the bins for the wavelet transform
+    Eq 4-6 in Wavelets paper
+    """
+    T = duration
+    N = data_len
+    fs = N / T
+    fmax = fs / 2
+
+    delta_t = T / Nt
+    delta_f = 1 / (2 * delta_t)
+
+    # assert delta_f == fmax / Nf, f"delta_f={delta_f} != fmax/Nf={fmax/Nf}"
+
+    f_bins = np.arange(0, Nf) * delta_f
+    t_bins = np.arange(0, Nt) * delta_t
+
+    return t_bins, f_bins
diff --git a/examples/lisa_wdb_snr.pdf b/examples/lisa_wdb_snr.pdf
new file mode 100644
index 0000000..a800f70
Binary files /dev/null and b/examples/lisa_wdb_snr.pdf differ
diff --git a/examples/lisa_wdb_snr.py b/examples/lisa_wdb_snr.py
new file mode 100644
index 0000000..d4b799c
--- /dev/null
+++ b/examples/lisa_wdb_snr.py
@@ -0,0 +1,162 @@
+from typing import Tuple
+
+import matplotlib.pyplot as plt
+
+# import TwoSlopeNorm from matplotlib.colors
+from matplotlib.colors import TwoSlopeNorm
+
+import numpy as np
+from scipy.signal.windows import tukey
+from common import (
+    compute_frequency_optimal_snr,
+    evolutionary_psd_from_stationary_psd,
+    compute_wavelet_snr,
+    get_wavelet_bins,
+)
+from WDMWaveletTransforms.wavelet_transforms import transform_wavelet_time
+
+
+from collections import namedtuple
+
+TIMESERIES = namedtuple("TimeSeries", ["data", "time"])
+
+
+def lisa_psd_func(f):
+    """
+    PSD obtained from:
+    Robson et al 2018, "LISA Sensitivity Curves"
+    https://arxiv.org/pdf/1803.01944.pdf
+
+    Removed galactic confusion noise. Non stationary effect.
+
+    The power spectrum -- not the TDI
+
+    """
+
+    L = 2.5 * 10**9  # Length of LISA arm
+    f0 = 19.09 * 10**-3
+
+    # Eq 10
+    Poms = ((1.5 * 10**-11) ** 2) * (
+        1 + ((2 * 10**-3) / f) ** 4
+    )  # Optical Metrology Sensor
+
+    # Eq 11
+    Pacc = (
+        (3 * 10**-15) ** 2
+        * (1 + (4 * 10**-3 / (10 * f)) ** 2)
+        * (1 + (f / (8 * 10**-3)) ** 4)
+    )  # Acceleration Noise
+
+    # Eq 13
+    PSD = (
+        (10 / (3 * L**2))
+        * (Poms + (4 * Pacc) / ((2 * np.pi * f)) ** 4)
+        * (1 + 0.6 * (f / f0) ** 2)
+    )  # PSD
+
+    PSD = np.ones(len(PSD)) * 1e-40
+    return PSD
+
+
+def waveform(a: float, f: float, fdot: float, t: np.ndarray, eps=0):
+    """
+    This is a function. It takes in a value of the amplitude $a$, frequency $f$ and frequency derivative $\dot{f}
+    and a time vector $t$ and spits out whatever is in the return function. Modify amplitude to improve SNR.
+    Modify frequency range to also affect SNR but also to see if frequencies of the signal are important
+    for the windowing method. We aim to estimate the parameters $a$, $f$ and $\dot{f}$.
+
+    h = a * sin(2 * pi * (ft + 0.5 * fdot * t^2))
+    Quadratic chirp signal.
+
+    """
+
+    return a * (np.sin((2 * np.pi) * (f * t + 0.5 * fdot * t**2)))
+
+
+def get_lisa_data():
+    """
+    This function is used to generate the data for the LISA detector. We use the waveform function to generate
+    the signal and then use the freq_PSD function to generate the PSD. We then use the FFT function to generate
+    the frequency domain waveform. We then compute the optimal SNR.
+    """
+
+    a_true = 5e-21
+    f_true = 1e-3
+    fdot_true = 1e-8
+
+    fs = 2 * f_true  # Sampling rate
+    delta_t = np.floor(0.01 / fs)  # Sampling interval -- largely oversampling here.
+    tmax = 120 * 60 * 60  # 120 hours
+    t = np.arange(0, tmax, delta_t)
+    ND = int(
+        2 ** (np.ceil(np.log2(len(t))))
+    )  # Round length of time series to a power of two.
+    t = np.arange(0, ND) * delta_t
+
+    h_signal_t = waveform(a_true, f_true, fdot_true, t)
+    freq = np.fft.fftfreq(ND, delta_t)[: ND // 2]
+    psd_vals = lisa_psd_func(freq)
+    h_signal_f = np.fft.fft(h_signal_t)[: ND // 2]
+    duration = delta_t * ND
+
+    # skip first element to avoid division by zero
+    freq = freq[1:]
+    psd_vals = psd_vals[1:]
+    h_signal_f = h_signal_f[1:]
+
+    snr = compute_frequency_optimal_snr(h_signal_f, psd_vals, duration)
+    return freq, t, h_signal_t, h_signal_f, psd_vals, snr
+
+
+def main():
+    np.random.seed(1234)
+
+    freq, t, h_signal_t, h_signal_f, psd, snr = get_lisa_data()
+    ND = len(t)
+    Nf = 256
+    Nt = ND // Nf
+    dt = t[1]
+    duration = ND * dt
+
+    h_time = TIMESERIES(h_signal_t, t)
+    h_wavelet = transform_wavelet_time(h_time.data, Nf=Nf, Nt=Nt) * dt * np.sqrt(2)
+    time_grid, freq_grid = get_wavelet_bins(duration, ND, Nf, Nt)
+    hf = h_signal_f.data
+
+    # h_wavelet = from_time_to_wavelet(h_time, Nt=Nt)
+    psd_wavelet = evolutionary_psd_from_stationary_psd(psd, freq, freq_grid, time_grid) * dt
+    snr2 = (4*dt/ND) * np.sum(np.conj(hf)*hf/psd)
+    wavelet_snr2 = compute_wavelet_snr(h_wavelet, psd_wavelet)**2
+
+    print(f"SNR2 in time domain (ND:{ND}): {snr2:.2f}")
+    print(f"SNR2 in wavelet domain (Nf{Nf}xNt{Nt}): {wavelet_snr2:.2f}")
+
+    fig, axes = plt.subplots(2, 1, figsize=(5, 8))
+    axes[0].loglog(freq, np.abs(h_signal_f), label="Signal")
+    axes[0].loglog(freq, psd, label="PSD")
+    axes[0].legend()
+    axes[0].set_xlabel("Frequency [Hz]")
+    axes[0].set_ylabel("PSD")
+    axes[0].set_title("Frequency domain signal")
+    axes[1].pcolor(
+        time_grid, freq_grid, h_wavelet.T, cmap="RdBu", norm=TwoSlopeNorm(vcenter=0)
+    )
+    axes[1].set_ylim(1e-4, 1e-2)
+    axes[1].set_xlabel("Time [s]")
+    axes[1].set_ylabel("Frequency [Hz]")
+    axes[0].text(
+        0.1, 0.85, f"SNR: {snr:.2f}", transform=axes[0].transAxes, fontsize="x-large"
+    )
+    axes[1].text(
+        0.1,
+        0.85,
+        f"SNR: {wavelet_snr2:.2f}",
+        transform=axes[1].transAxes,
+        fontsize="x-large",
+    )
+    fig.savefig("lisa_wdb_snr.pdf", dpi=300)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/lvk_cbc_snr.pdf b/examples/lvk_cbc_snr.pdf
new file mode 100644
index 0000000..1759c84
Binary files /dev/null and b/examples/lvk_cbc_snr.pdf differ
diff --git a/examples/lvk_cbc_snr.py b/examples/lvk_cbc_snr.py
new file mode 100644
index 0000000..e601b7d
--- /dev/null
+++ b/examples/lvk_cbc_snr.py
@@ -0,0 +1,213 @@
+""" Example script to
+1. Generate a CBC signal with a given SNR in LVK 04 noise
+2. Perform a wavelet transform on the signal
+3. Compute the SNR of the wavelet transform
+4. Plot the wavelet transform + time domain signal
+
+
+REQUIRES:
+- WDMWaveletTransforms, bilby[gw], numba
+"""
+
+from WDMWaveletTransforms.wavelet_transforms import transform_wavelet_time
+import bilby
+from typing import Tuple
+import numpy as np
+from common import (
+    compute_wavelet_snr,
+    evolutionary_psd_from_stationary_psd,
+    get_wavelet_bins,
+    compute_frequency_optimal_snr,
+)
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+
+from collections import namedtuple
+
+TIMESERIES = namedtuple("TimeSeries", ["data", "time"])
+FREQSERIES = namedtuple("FreqSeries", ["data", "freq"])
+
+from scipy.interpolate import interp1d
+
+DURATION = 4
+SAMPLING_FREQUENCY = 16384
+DT = 1 / SAMPLING_FREQUENCY
+MINIMUM_FREQUENCY = 20
+MAXIMUM_FREQUENCY = 1024
+
+CBC_GENERATOR = bilby.gw.WaveformGenerator(
+    duration=DURATION,
+    sampling_frequency=SAMPLING_FREQUENCY,
+    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+    waveform_arguments=dict(
+        waveform_approximant="IMRPhenomD",
+        reference_frequency=20.0,
+        minimum_frequency=MINIMUM_FREQUENCY,
+        maximum_frequency=MAXIMUM_FREQUENCY,
+    ),
+)
+
+GW_PARMS = dict(
+    mass_1=30,
+    mass_2=30,  # 2 mass parameters
+    a_1=0.1,
+    a_2=0.1,
+    tilt_1=0.0,
+    tilt_2=0.0,
+    phi_12=0.0,
+    phi_jl=0.0,  # 6 spin parameters
+    ra=1.375,
+    dec=-1.2108,
+    luminosity_distance=2000.0,
+    theta_jn=0.0,  # 7 extrinsic parameters
+    psi=2.659,
+    phase=1.3,
+    geocent_time=0,
+)
+
+
+def _get_ifo(t0=0.0, noise=True):
+    ifos = bilby.gw.detector.InterferometerList(["H1"])  # design sensitivity
+    if noise:
+        ifos.set_strain_data_from_power_spectral_densities(
+            sampling_frequency=SAMPLING_FREQUENCY,
+            duration=DURATION,
+            start_time=t0,
+        )
+    else:
+        ifos.set_strain_data_from_zero_noise(
+            sampling_frequency=SAMPLING_FREQUENCY,
+            duration=DURATION,
+            start_time=t0,
+        )
+    return ifos
+
+
+def inject_signal_in_noise(
+    mc, q=1, distance=1000.0, noise=True
+) -> Tuple[TIMESERIES, float]:
+    injection_parameters = GW_PARMS.copy()
+    (
+        injection_parameters["mass_1"],
+        injection_parameters["mass_2"],
+    ) = bilby.gw.conversion.chirp_mass_and_mass_ratio_to_component_masses(mc, q)
+    injection_parameters["luminosity_distance"] = distance
+    ifos = _get_ifo(injection_parameters["geocent_time"] + 1.5, noise=noise)
+    ifos.inject_signal(
+        waveform_generator=CBC_GENERATOR, parameters=injection_parameters
+    )
+    ifo: bilby.gw.detector.Interferometer = ifos[0]
+    snr = ifo.meta_data["matched_filter_SNR"]
+    ifo.time_array = np.linspace(0, DURATION, int(DURATION * SAMPLING_FREQUENCY))
+    timeseries = TIMESERIES(ifo.strain_data.time_domain_strain, ifo.time_array)
+    fmask = ifo.frequency_mask
+    f = ifo.strain_data.frequency_array[fmask]
+    hf = ifo.strain_data.frequency_domain_strain[fmask]
+    psd = ifo.power_spectral_density_array[fmask]
+    freqseries = FREQSERIES(hf, f)
+    snr = compute_frequency_optimal_snr(hf, psd, DURATION)
+    return timeseries, freqseries, np.abs(snr), psd
+
+
+def plot(
+    signal_t,
+    signal_f,
+    optimal_snr_ht,
+    psd,
+    signal_wavelet,
+    time_grid,
+    freq_grid,
+    psd_wavelet,
+    snr_wavelet,
+    DURATION,
+    MINIMUM_FREQUENCY,
+):
+    # plot Timeseries, wavelet signal and PSD
+    fig, ax = plt.subplots(2, 2, figsize=(12, 10))
+    ax[0, 0].plot(signal_t.time, signal_t.data)
+    ax[0, 0].set_xlim(0, DURATION)
+    ax[0, 0].set_xlabel("Time [s]")
+    ax[0, 0].set_ylabel("Strain")
+    ax[0, 0].set_title("Time domain signal")
+    # ANNOTATE SNR
+    ax[0, 0].text(
+        0.1,
+        0.85,
+        f"SNR: {optimal_snr_ht:.2f}",
+        transform=ax[0, 0].transAxes,
+        fontsize="x-large",
+    )
+    ax[0, 1].loglog(signal_f.freq, psd)
+    ax[0, 1].loglog(signal_f.freq, np.abs(signal_f.data))
+
+    ax[0, 1].set_title("PSD in frequency domain")
+
+    ax[1, 0].pcolor(
+        time_grid,
+        freq_grid,
+        signal_wavelet.T,
+        cmap="RdBu",
+        norm=colors.TwoSlopeNorm(vcenter=0),
+    )
+    ax[1, 0].set_ylim(MINIMUM_FREQUENCY, MAXIMUM_FREQUENCY / 2)
+    ax[1, 0].set_xlabel("Time [s]")
+    ax[1, 0].set_ylabel("Frequency [Hz]")
+    ax[1, 0].set_title("Wavelet transform")
+    ax[1, 0].text(
+        0.1,
+        0.85,
+        f"SNR: {snr_wavelet:.2f}",
+        transform=ax[1, 0].transAxes,
+        fontsize="x-large",
+    )
+    ax[1, 1].pcolor(time_grid, freq_grid, np.log(psd_wavelet.T), cmap="viridis")
+    ax[1, 1].set_ylim(MINIMUM_FREQUENCY, MAXIMUM_FREQUENCY / 2)
+    ax[1, 1].set_title("PSD in wavelet domain")
+    ax[1, 1].set_xlabel("Time [s]")
+    ax[1, 1].set_ylabel("Frequency [Hz]")
+    return fig
+
+
+def main():
+    signal_t, signal_f, optimal_snr_ht, psd = inject_signal_in_noise(mc=30, noise=False)
+    ND = len(signal_t.time)
+    Nt = 128
+    Nf = ND // Nt
+    dt = signal_t.time[1] - signal_t.time[0]
+
+    signal_wavelet = transform_wavelet_time(signal_t.data, Nf, Nt) * dt * np.sqrt(2)
+    time_grid, freq_grid = get_wavelet_bins(DURATION, ND, Nf, Nt)
+    psd_wavelet = (
+        evolutionary_psd_from_stationary_psd(psd, signal_f.freq, freq_grid, time_grid)
+        * dt
+    )
+
+    # ignore the freq below 20 Hz
+    # get idx where freq_grid > 20hz
+    idx = np.where(freq_grid > MINIMUM_FREQUENCY)[0][0]
+    freq_grid = freq_grid[idx:]
+    psd_wavelet = psd_wavelet[:, idx:]
+    signal_wavelet = signal_wavelet[:, idx:]
+    snr_wavelet = compute_wavelet_snr(signal_wavelet, psd_wavelet)
+    print("SNR in time domain:", optimal_snr_ht)
+    print(f"ND: {ND}, Nt: {Nt}, Nf: {Nf}")
+    print("SNR in wavelet domain:", snr_wavelet)
+    plot(
+        signal_t,
+        signal_f,
+        optimal_snr_ht,
+        psd,
+        signal_wavelet,
+        time_grid,
+        freq_grid,
+        psd_wavelet,
+        snr_wavelet,
+        DURATION,
+        MINIMUM_FREQUENCY,
+    ).savefig("lvk_cbc_snr.pdf", dpi=300)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/toy_model_snr.pdf b/examples/toy_model_snr.pdf
new file mode 100644
index 0000000..85efc9c
Binary files /dev/null and b/examples/toy_model_snr.pdf differ
diff --git a/examples/toy_model_snr.py b/examples/toy_model_snr.py
new file mode 100644
index 0000000..3952330
--- /dev/null
+++ b/examples/toy_model_snr.py
@@ -0,0 +1,146 @@
+"""
+
+Toy model of sine-wave signal and a constant PSD
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+from WDMWaveletTransforms.wavelet_transforms import transform_wavelet_time, inverse_wavelet_time
+import pytest
+import scipy
+
+def flat_psd_func(f, psd_amp):
+    return psd_amp * np.ones(len(f))
+
+
+def colored_psd_func(f, psd_amp):
+    return psd_amp * np.ones(len(f)) * (1 + (f / 10) ** 2)
+
+
+def colored_psd_func2(f, psd_amp):
+    return psd_amp * np.ones(len(f)) * (1 + 1 / (f / 10) ** 2)
+
+
+# pytest parameterize decorator
+@pytest.mark.parametrize(
+    "f0, T, A, PSD_AMP, Nf, psd_func",
+    [
+        (20, 1000, 1e-3, 1e-2, 16, flat_psd_func),
+        (10, 1000, 1e-3, 1e-2, 32, colored_psd_func2),
+        (20, 1000, 1e-3, 1e-2, 16, colored_psd_func),
+    ])
+def test_wavelet_timedomain_snr(f0, T, A, PSD_AMP, Nf, psd_func):
+    ########################################
+    # Part1: Analytical SNR calculation
+    ########################################
+    dt = 0.5 / (2 * f0)  # Shannon's sampling theorem, set dt < 1/2*highest_freq
+    t = np.arange(0, T, dt)  # Time array
+    # round len(t) to the nearest power of 2
+    t = t[:2 ** int(np.log2(len(t)))]
+    T = len(t) * dt
+
+    y = A * np.sin(2 * np.pi * f0 * t)  # Signal waveform we wish to test
+
+    freq = np.fft.fftshift(np.fft.fftfreq(len(y), dt))  # Frequencies
+    df = abs(freq[1] - freq[0])  # Sample spacing in frequency
+
+    y_fft = dt * np.fft.fftshift(np.fft.fft(y))  # continuous time fourier transform [seconds]
+    N_f = len(y_fft)
+    N_t = len(y)
+
+    PSD = psd_func(freq, PSD_AMP)  # PSD of the noise
+
+    # Compute the SNRs
+    SNR2_f = 2 * np.sum(abs(y_fft) ** 2 / PSD) * df
+    SNR2_t = 2 * dt * np.sum(abs(y) ** 2 / PSD)
+    SNR2_t_analytical = (A ** 2) * T / PSD[0]
+
+    ########################################
+    # Part2: Wavelet domain
+    ########################################
+
+    ND = len(y)
+    Nt = ND // Nf
+    ND = Nf * Nt
+
+    signal_wavelet = transform_wavelet_time(y, Nf=Nf, Nt=Nt) * np.sqrt(2) * dt
+
+    delta_t = T / Nt
+    delta_f = 1 / (2 * delta_t)
+    freq_grid = np.arange(0, Nf) * delta_f
+    time_grid = np.arange(0, Nt) * delta_t
+    psd = psd_func(freq_grid, PSD_AMP)
+    psd_wavelet = np.repeat(psd[None, :], Nt, axis=0) * dt
+
+    wavelet_snr2 = np.sum((signal_wavelet * signal_wavelet / psd_wavelet))
+    mse = np.mean((y - inverse_wavelet_time(signal_wavelet, Nf=Nf, Nt=Nt)) ** 2)
+    print('---------')
+    print(f"SNR squared in the frequency domain is = {SNR2_f:.2f}")
+    print(f"SNR squared in the time domain (Parseval's theorem) is = {SNR2_t:.2f}", )
+    print(f"(pen and paper) Analytical result would predict SNR squared = {SNR2_t_analytical:.2f}")
+    print(f"In the wavelet domain, SNR_sqr = {wavelet_snr2:.2f}")
+    print(f"Mean squared error in the wavelet domain = {mse:.2f}")
+    print('---------')
+    assert np.isclose(SNR2_f, wavelet_snr2, atol=1e-2), "SNR in time domain and wavelet domain should be the same"
+
+
+
+def test_chirp_signal():
+    T = 1000
+    f0 = 1e-3
+    f1 = 20
+    dt = 0.5 / (2 * f1)  # Shannon's sampling theorem, set dt < 1/2*highest_freq
+    t = np.arange(0, T, dt)  # Time array
+    # round len(t) to the nearest power of 2
+    t = t[:2 ** int(np.log2(len(t)))]
+    T = len(t) * dt
+    y = scipy.signal.chirp(t, f0=f0, f1=f1, t1=T, method='quadratic')
+    # plot spectogram
+    plt.specgram(y, Fs=1/dt, NFFT=128, noverlap=64, cmap='viridis')
+    plt.show()
+
+
+    freq = np.fft.fftshift(np.fft.fftfreq(len(y), dt))  # Frequencies
+    df = abs(freq[1] - freq[0])  # Sample spacing in frequency
+
+    y_fft = dt * np.fft.fftshift(np.fft.fft(y))  # continuous time fourier transform [seconds]
+    N_f = len(y_fft)
+    N_t = len(y)
+
+    psd_func = flat_psd_func
+    PSD_AMP = 1e-2
+    PSD = psd_func(freq, PSD_AMP)  # PSD of the noise
+
+    # Compute the SNRs
+    SNR2_f = 2 * np.sum(abs(y_fft) ** 2 / PSD) * df
+
+    ########################################
+    # Part2: Wavelet domain
+    ########################################
+    Nf = 16
+    ND = len(y)
+    Nt = ND // Nf
+    ND = Nf * Nt
+
+    signal_wavelet = transform_wavelet_time(y, Nf=Nf, Nt=Nt) * np.sqrt(2) * dt
+
+    delta_t = T / Nt
+    delta_f = 1 / (2 * delta_t)
+    freq_grid = np.arange(0, Nf) * delta_f
+    time_grid = np.arange(0, Nt) * delta_t
+    psd = psd_func(freq_grid, PSD_AMP)
+    psd_wavelet = np.repeat(psd[None, :], Nt, axis=0) * dt
+
+    wavelet_snr2 = np.sum((signal_wavelet * signal_wavelet / psd_wavelet))
+    mse = np.mean((y - inverse_wavelet_time(signal_wavelet, Nf=Nf, Nt=Nt)) ** 2)
+    print('---------')
+
+    assert np.isclose(SNR2_f, wavelet_snr2, atol=1e-2), "SNR in time domain and wavelet domain should be the same"
+
+
+
+
+
+
+if __name__ == '__main__':
+    test_chirp_signal()