diff --git a/pyadjoint/adjoint_source_types/ccc_misfit.py b/pyadjoint/adjoint_source_types/ccc_misfit.py
new file mode 100644
index 0000000..a09d3d6
--- /dev/null
+++ b/pyadjoint/adjoint_source_types/ccc_misfit.py
@@ -0,0 +1,123 @@
+"""
+Cross correlation coefficient waveform misfit and adjoint source.
+
+This file will also serve as an explanation of how to add new adjoint
+sources to Pyadjoint.
+
+:copyright:
+    Lion Krischer (krischer@geophysik.uni-muenchen.de), 2015
+    Shi Yao (yaoshi229@gmail.com)
+:license:
+    BSD 3-Clause ("BSD New" or "BSD Simplified")
+"""
+import numpy as np
+from pyadjoint.utils.signal import window_taper
+
+def calculate_adjoint_source(observed, synthetic, config, windows,
+                             observed_2=None, synthetic_2=None, windows_2=None, 
+                             adjoint_src=True, window_stats=True):
+    """
+    Calculate adjoint source for the Cross correlation coefficient misfit
+    measurement
+
+    :type observed: obspy.core.trace.Trace
+    :param observed: observed waveform to calculate adjoint source
+    :type synthetic:  obspy.core.trace.Trace
+    :param synthetic: synthetic waveform to calculate adjoint source
+    :type config: pyadjoint.config.ConfigCCC
+    :param config: Config class with parameters to control processing
+    :type windows: list of tuples
+    :param windows: [(left, right),...] representing left and right window
+        borders to be used to calculate misfit and adjoint sources
+    :type observed_2: obspy.core.trace.Trace
+    :param observed_2: second observed waveform to calculate adjoint sources
+        from station pairs
+    :type synthetic_2:  obspy.core.trace.Trace
+    :param synthetic_2: second synthetic waveform to calculate adjoint sources
+        from station pairs
+    :type windows_2: list of tuples
+    :param windows_2: [(left, right),...] representing left and right window
+        borders to be tapered in units of seconds since first sample in data
+        array. Used to window `observed_2` and `synthetic_2`
+    """
+    assert(config.__class__.__name__ == "ConfigTFPhase"), \
+        "Incorrect configuration class passed to Time-Frequency Phase misfit"
+
+    if config.double_difference:
+        raise NotImplementedError(
+            "Time-Frequency phase misfit does not have double difference "
+            "capabilities"
+        )
+    
+    # Dictionary of return values related to exponentiated phase
+    ret_val = {}
+
+    # List of windows and some measurement values for each
+    win_stats = []
+
+    nlen_data = len(synthetic.data)
+    deltat = synthetic.stats.delta
+
+    adj = np.zeros(nlen_data)
+
+    misfit_sum = 0.0
+
+    # loop over time windows
+    for window in windows:
+
+        left_window_border = window[0]
+        right_window_border = window[1]
+
+        left_sample = int(np.floor(left_window_border / deltat)) + 1
+        nlen = int(np.floor((right_window_border - left_window_border) /
+                            deltat)) + 1
+        right_sample = left_sample + nlen
+
+        d = np.zeros(nlen)
+        s = np.zeros(nlen)
+
+        d[0: nlen] = observed.data[left_sample: right_sample]
+        s[0: nlen] = synthetic.data[left_sample: right_sample]
+
+        # All adjoint sources will need some kind of windowing taper
+        # to get rid of kinks at two ends
+        window_taper(d, taper_percentage=config.taper_percentage,
+                     taper_type=config.taper_type)
+        window_taper(s, taper_percentage=config.taper_percentage,
+                     taper_type=config.taper_type)
+
+        #diff = s - d
+        CC = np.dot(d,s)
+        weight_2 = np.sqrt(np.dot(d, d) * np.dot(s, s))
+        misfit = 1 - CC / weight_2
+        
+        A = np.dot(d, s) / np.dot(s, s)
+        # print(A)
+        diff_w = (d - A * s) / weight_2
+        #diff_w = diff * -1.0
+        window_taper(diff_w, taper_percentage=config.taper_percentage,
+                     taper_type=config.taper_type)
+        # for some reason the 0.5 (see 2012 measure_adj mannual, P11) is
+        # not in misfit definetion in measure_adj
+        # misfit_sum += 0.5 * simps(y=diff_w**2, dx=deltat)
+        misfit_sum += misfit
+
+        adj[left_sample: right_sample] = diff_w
+
+        win_stats.append(
+            {"type": config.adjsrc_type, "measurement_type": "dt",
+             "left": left_sample * deltat, "right": right_sample * deltat,
+             "misfit" : misfit
+            }
+        )
+    
+    ret_val["misfit"] = misfit_sum
+
+    if window_stats:
+        ret_val["window_stats"] = win_stats
+
+    if adjoint_src:
+        # Reverse in time
+        ret_val["adjoint_source"] = -adj[::-1]
+
+    return ret_val        
diff --git a/pyadjoint/adjoint_source_types/tf_phase_misfit.py b/pyadjoint/adjoint_source_types/tf_phase_misfit.py
new file mode 100644
index 0000000..bffac4f
--- /dev/null
+++ b/pyadjoint/adjoint_source_types/tf_phase_misfit.py
@@ -0,0 +1,289 @@
+"""
+An implementation of the time frequency phase misfit and adjoint source after
+Fichtner et al. (2008). This is different from lasif version in order to account for
+structure of the recent version of pyadjoint.
+
+:copyright:
+    Lion Krischer (krischer@geophysik.uni-muenchen.de), 2013
+    Yajian Gao (krischer@geophysik.uni-muenchen.de), 2021
+    Shi Yao (yaoshi229@gmail.com), 2024
+:license:
+    GNU General Public License, Version 3
+    (http://www.gnu.org/copyleft/gpl.html)
+"""
+import warnings
+
+import numexpr as ne
+import numpy as np
+from obspy.signal.interpolation import lanczos_interpolation
+from pyadjoint.utils.signal import get_window_info, window_taper
+from pyadjoint.timefrequency_utils import utils
+from pyadjoint.timefrequency_utils import time_frequency
+
+
+def calculate_adjoint_source(observed, synthetic, config, windows,
+                             observed_2=None, synthetic_2=None, windows_2=None, 
+                             adjoint_src=True, window_stats=True):
+    """
+    Calculate adjoint source for the time-frequency phase misfit
+    measurement
+
+    :type observed: obspy.core.trace.Trace
+    :param observed: observed waveform to calculate adjoint source
+    :type synthetic:  obspy.core.trace.Trace
+    :param synthetic: synthetic waveform to calculate adjoint source
+    :type config: pyadjoint.config.ConfigTFPhase
+    :param config: Config class with parameters to control processing
+    :type windows: list of tuples
+    :param windows: [(left, right),...] representing left and right window
+        borders to be used to calculate misfit and adjoint sources
+    :type observed_2: obspy.core.trace.Trace
+    :param observed_2: second observed waveform to calculate adjoint sources
+        from station pairs
+    :type synthetic_2:  obspy.core.trace.Trace
+    :param synthetic_2: second synthetic waveform to calculate adjoint sources
+        from station pairs
+    :type windows_2: list of tuples
+    :param windows_2: [(left, right),...] representing left and right window
+        borders to be tapered in units of seconds since first sample in data
+        array. Used to window `observed_2` and `synthetic_2`
+    """
+    assert(config.__class__.__name__ == "ConfigTFPhase"), \
+        "Incorrect configuration class passed to Time-Frequency Phase misfit"
+
+    if config.double_difference:
+        raise NotImplementedError(
+            "Time-Frequency phase misfit does not have double difference "
+            "capabilities"
+        )
+
+    # Dictionary of return values related to exponentiated phase
+    ret_val = {}
+
+    # List of windows and some measurement values for each
+    win_stats = []
+
+    # Assumes that t starts at 0. Pad your data if that is not the case -
+    # Parts with zeros are essentially skipped making it fairly efficient.
+    t = observed.times(type="relative")
+    max_criterion = 7.0
+    taper = True
+    taper_ratio = 0.3
+    taper_type="cosine"
+    assert t[0] == 0
+
+    window_weight = 1.0
+
+    # Initiate constants and empty return values to fill
+    nlen_w_data = len(synthetic.data)
+    dt = synthetic.stats.delta
+
+    adj = np.zeros(nlen_w_data)  # adjoint source
+
+    misfit_sum = 0.0
+
+    # loop over time windows
+    for window in windows:
+        left_sample, right_sample, nlen_w = get_window_info(window, dt)
+
+        # Initiate empty window arrays for memory efficiency
+        observed = utils.window_trace(
+            trace=observed,
+            window=window,
+            taper=taper,
+            taper_ratio=taper_ratio,
+            taper_type=taper_type
+        )
+        synthetic = utils.window_trace(
+            trace=synthetic,
+            window=window,
+            taper=taper,
+            taper_ratio=taper_ratio,
+            taper_type=taper_type
+        )
+        
+        # Internal sampling interval. Some explanations for this "magic" number.
+        # LASIF's preprocessing allows no frequency content with smaller periods
+        # than min_period / 2.2 (see function_templates/preprocesssing_function.py
+        # for details). Assuming most users don't change this, this is equal to
+        # the Nyquist frequency and the largest possible sampling interval to
+        # catch everything is min_period / 4.4.
+        #
+        # The current choice is historic as changing does (very slightly) chance
+        # the calculated misfit and we don't want to disturb inversions in
+        # progress. The difference is likely minimal in any case. We might have
+        # same aliasing into the lower frequencies but the filters coupled with
+        # the TF-domain weighting will get rid of them in essentially all
+        # realistically occurring cases.
+        dt_new = max(float(int(config.min_period / 4.0)), t[1] - t[0])
+        dt_old = t[1] - t[0]
+
+        # new time axis
+        ti = utils.matlab_range(t[0], t[-1], dt_new)
+        # Make sure its odd - that avoid having to deal with some issues
+        # regarding frequency bin interpolation. Now positive and negative
+        # frequencies will always be all symmetric. Data is assumed to be
+        # tapered in any case so no problem are to be expected.
+        if not len(ti) % 2:
+            ti = ti[:-1]
+
+        # Interpolate both signals to the new time axis - this massively speeds
+        # up the whole procedure as most signals are highly oversampled. The
+        # adjoint source at the end is re-interpolated to the original sampling
+        # points.        
+        data = lanczos_interpolation(
+            data=observed.data,
+            old_start=t[0],
+            old_dt=t[1] - t[0],
+            new_start=t[0],
+            new_dt=dt_new,
+            new_npts=len(ti),
+            a=8,
+            window="blackmann",
+            )
+        synthetic_data = lanczos_interpolation(
+            data=synthetic.data,
+            old_start=t[0],
+            old_dt=t[1] - t[0],
+            new_start=t[0],
+            new_dt=dt_new,
+            new_npts=len(ti),
+            a=8,
+            window="blackmann",
+            )
+        original_time = t
+        t = ti        
+
+        # -------------------------------------------------------------------------
+        # Compute time-frequency representations
+
+        # Window width is twice the minimal period.
+        width = 2.0 * config.min_period
+
+        # Compute time-frequency representation of the cross-correlation
+        _, _, tf_cc = time_frequency.time_frequency_cc_difference(
+            t, data, synthetic_data, width
+            )
+        # Compute the time-frequency representation of the synthetic
+        tau, nu, tf_synth = time_frequency.time_frequency_transform(
+            t, synthetic_data, width
+            )
+        
+        # -------------------------------------------------------------------------
+        # compute tf window and weighting function
+
+        # noise taper: down-weight tf amplitudes that are very low
+        tf_cc_abs = np.abs(tf_cc)
+        m = tf_cc_abs.max() / 10.0  # NOQA
+        weight = ne.evaluate("1.0 - exp(-(tf_cc_abs ** 2) / (m ** 2))")
+        nu_t = nu.T 
+
+        # highpass filter (periods longer than max_period are suppressed
+        # exponentially)
+        weight *= 1.0 - np.exp(-((nu_t * config.max_period) ** 2))  
+
+        # lowpass filter (periods shorter than min_period are suppressed
+        # exponentially)
+        nu_t_large = np.zeros(nu_t.shape)
+        nu_t_small = np.zeros(nu_t.shape)
+        thres = nu_t <= 1.0 / config.min_period
+        nu_t_large[np.invert(thres)] = 1.0
+        nu_t_small[thres] = 1.0
+        weight *= (
+            np.exp(-10.0 * np.abs(nu_t * config.min_period - 1.0)) * nu_t_large
+            + nu_t_small
+            )
+
+        # normalization
+        weight /= weight.max()                     
+
+        # computation of phase difference, make quality checks and misfit ---------
+
+        # Compute the phase difference.
+        # DP = np.imag(np.log(m + tf_cc / (2 * m + np.abs(tf_cc))))
+        DP = np.angle(tf_cc)
+
+        # Attempt to detect phase jumps by taking the derivatives in time and
+        # frequency direction. 0.7 is an emperical value.
+        abs_weighted_DP = np.abs(weight * DP)
+        _x = abs_weighted_DP.max()  # NOQA
+        test_field = ne.evaluate("weight * DP / _x")
+
+        criterion_1 = np.sum([np.abs(np.diff(test_field, axis=0)) > 0.7])
+        criterion_2 = np.sum([np.abs(np.diff(test_field, axis=1)) > 0.7])
+        criterion = np.sum([criterion_1, criterion_2])
+
+        # Compute the phase misfit
+        dnu = nu[1] - nu[0]          
+
+        i = ne.evaluate("sum(weight ** 2 * DP ** 2)")
+
+        phase_misfit = np.sqrt(i * dt_new * dnu) * window_weight
+
+        misfit_sum += phase_misfit
+
+        if np.isnan(phase_misfit):
+            print("The phase misfit is NaN.")
+            raise ValueError("Phase misfit cannot be NaN.")
+        
+        # The misfit can still be computed, even if not adjoint source is
+        # available.
+        if criterion > max_criterion:
+            warning = (
+            "Possible phase jump detected. Misfit included. No "
+            "adjoint source computed. Criterion: %.1f - Max allowed "
+            "criterion: %.1f" % (criterion, max_criterion)
+            )
+            warnings.warn(warning)
+
+        # Make kernel for the inverse tf transform
+        idp = ne.evaluate(
+            "weight ** 2 * DP * tf_synth / (m + abs(tf_synth) ** 2)"
+            )
+        
+        # Invert tf transform and make adjoint source
+        ad_src, it, I = time_frequency.itfa(tau, idp, width)
+
+        # Interpolate both signals to the new time axis
+        # Pad with a couple of zeros in case some where lost in all
+        # these resampling operations. The first sample should not
+        # change the time.
+        ad_src = lanczos_interpolation(
+                data=np.concatenate([ad_src.imag, np.zeros(100)]),
+                old_start=tau[0],
+                old_dt=tau[1] - tau[0],
+                new_start=original_time[0],
+                new_dt=original_time[1] - original_time[0],
+                new_npts=len(original_time),
+                a=8,
+                window="blackmann",
+            )        
+                     
+        # Divide by the misfit and change sign.
+        ad_src /= phase_misfit + np.spacing(1)
+        ad_src = ad_src / ((t[1] - t[0]) ** 2) * dt_old
+        
+        # taper the adjoint source
+        ad_src_short = ad_src[left_sample: right_sample]
+        window_taper(ad_src_short, taper_percentage=config.taper_percentage,
+            taper_type=config.taper_type)
+        adj[left_sample: right_sample] = ad_src_short * window_weight        
+
+        win_stats.append(
+            {"type": config.adjsrc_type, "measurement_type": "dt",
+             "left": left_sample * dt, "right": right_sample * dt,
+             "misfit" : phase_misfit
+            }
+        )
+
+    ret_val["misfit"] = misfit_sum
+
+    if window_stats:
+        ret_val["window_stats"] = win_stats
+    
+    if adjoint_src is True:
+        # Reverse time and add a leading zero so the adjoint source has the
+        # same length as the input time series.
+        ret_val["adjoint_source"] = -adj[::-1]
+        
+    return ret_val
diff --git a/pyadjoint/config.py b/pyadjoint/config.py
index 9ea493a..67dc19a 100644
--- a/pyadjoint/config.py
+++ b/pyadjoint/config.py
@@ -18,10 +18,11 @@
 
 # Constants defining the available adjoint source types in this package and
 # their verbose names. New adjoint sources need to be added here.
+# REVISED by YAOSHI
 ADJSRC_TYPES = [
     "waveform", "convolution", "exponentiated_phase", "cc_traveltime",
-    "multitaper", "waveform_dd", "convolution_dd", "cc_traveltime_dd",
-    "multitaper_dd",
+    "multitaper", "tf_phase", "ccc", "waveform_dd", "convolution_dd", "cc_traveltime_dd",
+    "multitaper_dd", "tf_phase_dd", "ccc_dd"
     # >>> ADD NEW ADJOINT SOURCES HERE
 ]
 
@@ -55,6 +56,13 @@ def get_config(adjsrc_type, min_period, max_period, **kwargs):
     elif adjsrc_type in ["multitaper", "multitaper_dd"]:
         cfg = ConfigMultitaper(min_period, max_period, double_difference=dd,
                                **kwargs)
+    # REVISED by YAOSHI
+    elif adjsrc_type in ["tf_phase", "tf_phase_dd"]:
+        cfg = ConfigTFPhase(min_period, max_period, double_difference=dd,
+                               **kwargs)
+    elif adjsrc_type in ["ccc", "ccc_dd"]:
+        cfg = ConfigTFPhase(min_period, max_period, double_difference=dd,
+                                **kwargs)
     # >>> ADD NEW ADJOINT SOURCES HERE
     else:
         raise NotImplementedError(f"adjoint source type must be in "
@@ -105,6 +113,15 @@ def get_function(adjsrc_type):
         from pyadjoint.adjoint_source_types.multitaper_misfit import \
             calculate_adjoint_source
         fct = calculate_adjoint_source
+    # REVISED by YAOSHI
+    elif adjsrc_type in ["tf_phase", "tf_phase_dd"]:
+        from pyadjoint.adjoint_source_types.tf_phase_misfit import \
+            calculate_adjoint_source
+        fct = calculate_adjoint_source
+    elif adjsrc_type in ["ccc", "ccc_dd"]:
+        from pyadjoint.adjoint_source_types.ccc_misfit import \
+            calculate_adjoint_source
+        fct = calculate_adjoint_source
     # >>> ADD NEW ADJOINT SOURCES HERE
 
     return fct
@@ -317,4 +334,66 @@ def __init__(self, min_period, max_period, lnpt=15,
         # To be overwritten by get_config()
         self.adjsrc_type = None
 
+class ConfigTFPhase:
+    """
+    Time Frequency Phase misfit function required parameters
 
+    :param min_period: Minimum period of the filtered input data in seconds.
+    :type min_period: float
+    :param max_period: Maximum period of the filtered input data in seconds.
+    :type max_period: float
+    :param taper_percentage: Percentage of a time window needs to be
+        tapered at two ends, to remove the non-zero values for adjoint
+        source and for fft.
+    :type taper_percentage: float
+    :param taper_type: taper type, see pyadjoint.utils.signal.TAPER_COLLECTION
+        list for available taper types
+    :type taper_type: str
+    :param wtr_env: float
+    :param wtr_env: window taper envelope amplitude scaling
+    :type double_difference: bool
+    :param double_difference: flag to turn on double difference measurements,
+        which signals to the main calc function whether additional waveforms
+        are required at input
+    """
+    def __init__(self, min_period, max_period, taper_type="hann",
+                 taper_percentage=0.3, double_difference=False):
+        self.min_period = min_period
+        self.max_period = max_period
+        self.taper_type = taper_type
+        self.taper_percentage = taper_percentage
+        self.double_difference = double_difference
+        # To be overwritten by get_config()
+        self.adjsrc_type = None
+
+class ConfigCCC:
+    """
+    Cross correlation coefficient misfit function required parameters
+
+    :param min_period: Minimum period of the filtered input data in seconds.
+    :type min_period: float
+    :param max_period: Maximum period of the filtered input data in seconds.
+    :type max_period: float
+    :param taper_percentage: Percentage of a time window needs to be
+        tapered at two ends, to remove the non-zero values for adjoint
+        source and for fft.
+    :type taper_percentage: float
+    :param taper_type: taper type, see pyadjoint.utils.signal.TAPER_COLLECTION
+        list for available taper types
+    :type taper_type: str
+    :param wtr_env: float
+    :param wtr_env: window taper envelope amplitude scaling
+    :type double_difference: bool
+    :param double_difference: flag to turn on double difference measurements,
+        which signals to the main calc function whether additional waveforms
+        are required at input
+    """
+    def __init__(self, min_period, max_period, taper_type="hann",
+                 taper_percentage=0.3, double_difference=False):
+        self.min_period = min_period
+        self.max_period = max_period
+        self.taper_type = taper_type
+        self.taper_percentage = taper_percentage
+        self.double_difference = double_difference
+        # To be overwritten by get_config()
+        self.adjsrc_type = None
diff --git a/pyadjoint/main.py b/pyadjoint/main.py
index 47743b9..05336f3 100644
--- a/pyadjoint/main.py
+++ b/pyadjoint/main.py
@@ -221,10 +221,16 @@ def plot_adjoint_source(observed, synthetic, adjoint_source, misfit, windows,
     right_window_border = min(x_range, right_window_border)
 
     plt.subplot(211)
-    plt.plot(observed.times(), observed.data, color="0.2", label="Observed",
+    # REVISED by YAOSHI add the normalization for data
+    print(observed.times().min(), observed.times().max(), len(observed.times()))
+    plt.plot(observed.times(), observed.data / np.abs(observed.data).max(), color="0.2", label="Observed",
              lw=2)
-    plt.plot(synthetic.times(), synthetic.data, color="#bb474f",
+    plt.plot(synthetic.times(), synthetic.data / np.abs(synthetic.data).max(), color="#bb474f",
              label="Synthetic", lw=2)
+    # plt.plot(observed.times(), observed.data / np.abs(observed.data).max(), color="0.2", label="Observed",
+    #          lw=2)
+    # plt.plot(synthetic.times(), synthetic.data / np.abs(synthetic.data).max(), color="#bb474f",
+    #          label="Synthetic", lw=2)
     for window in windows:
         re = patches.Rectangle((window[0], plt.ylim()[0]),
                                window[1] - window[0],
@@ -234,25 +240,31 @@ def plot_adjoint_source(observed, synthetic, adjoint_source, misfit, windows,
 
     plt.grid()
     plt.legend(fancybox=True, framealpha=0.5)
-    plt.xlim(left_window_border, right_window_border)
+    # REVISED by YAOSHI
+    # plt.xlim(left_window_border, right_window_border)
 
+    # REVISED by YAOSHI
     # Determine min and max amplitudes within the time window
-    obs_win = observed.data[int(left_window_border):int(right_window_border)]
-    syn_win = synthetic.data[int(left_window_border):int(right_window_border)]
-    ylim = max([max(np.abs(obs_win)), max(np.abs(syn_win))])
-    plt.ylim(-ylim, ylim)
+    # obs_win = observed.data[int(left_window_border):int(right_window_border)]
+    # syn_win = synthetic.data[int(left_window_border):int(right_window_border)]
+    # ylim = max([max(np.abs(obs_win)), max(np.abs(syn_win))])
+    # plt.ylim(-ylim, ylim)
 
     plt.subplot(212)
-    plt.plot(observed.times(), adjoint_source[::-1], color="#2f8d5b", lw=2,
-             label="Adjoint Source")
+    # REVISED by YAOSHI, Delete the reverse, to keep the figure consisitent with the data
+    plt.plot(observed.times(), adjoint_source / np.abs(adjoint_source).max(), color="#2f8d5b", lw=2,
+            label="Adjoint Source")
+    # plt.plot(observed.times(), adjoint_source[::-1], color="#2f8d5b", lw=2,
+    #          label="Adjoint Source")
     plt.grid()
     plt.legend(fancybox=True, framealpha=0.5)
 
     # No time reversal for comparison with data
-    plt.xlim(left_window_border, right_window_border)
+    # Revised by YAOSHI, comment the xlim and ylim to show all data
+    # plt.xlim(left_window_border, right_window_border)
     plt.xlabel("Time (seconds)")
-    ylim = max(map(abs, plt.ylim()))
-    plt.ylim(-ylim, ylim)
+    # ylim = max(map(abs, plt.ylim()))
+    # plt.ylim(-ylim, ylim)
 
     plt.suptitle("%s Adjoint Source with a Misfit of %.3g" % (
         adjoint_source_name, misfit))
diff --git a/pyadjoint/timefrequency_utils/time_frequency.py b/pyadjoint/timefrequency_utils/time_frequency.py
new file mode 100644
index 0000000..12e0503
--- /dev/null
+++ b/pyadjoint/timefrequency_utils/time_frequency.py
@@ -0,0 +1,118 @@
+import numpy as np
+import scipy.fftpack
+import scipy.interpolate
+from pyadjoint.timefrequency_utils import utils
+
+
+def time_frequency_transform(t, s, width, threshold=1e-2):
+    """
+    Gabor transform (time frequency transform with Gaussian windows).
+
+    Data will be resampled before it is transformed.
+
+    :param t: discrete time.
+    :param s: discrete signal.
+    :param width: width of the Gaussian window
+    :param threshold: fraction of the absolute signal below which the Fourier
+        transform is set to zero in order to reduce computation time
+    """
+    N = len(t)
+    dt = t[1] - t[0]
+
+    nu = np.linspace(0, float(N - 1) / (N * dt), N)
+
+    # Compute the time frequency representation
+    tfs = np.zeros((N, N), dtype="complex128")
+
+    threshold = np.abs(s).max() * threshold
+
+    for k in range(N):
+        # Window the signals
+        f = utils.gaussian_window(t - t[k], width) * s
+
+        # No need to transform if nothing is there. Great speedup as lots of
+        # windowed functions have 0 everywhere.
+        if np.abs(f).max() < threshold:
+            continue
+
+        tfs[k, :] = scipy.fftpack.fft(f)
+
+    tfs *= dt / np.sqrt(2.0 * np.pi)
+
+    return t, nu, tfs
+
+
+def time_frequency_cc_difference(t, s1, s2, width, threshold=1e-2):
+    """
+    Straight port of tfa_cc_new.m
+
+    :param t: discrete time
+    :param s1: discrete signal 1
+    :param s2: discrete signal 2
+    :param width: width of the Gaussian window
+    :param threshold: fraction of the absolute signal below which the Fourier
+        transform is set to zero in order to reduce computation time
+    """
+    dt = t[1] - t[0]
+
+    # Extend the time axis, required for the correlation
+    N = len(t)
+    t_cc = np.linspace(t[0], t[0] + (2 * N - 2) * dt, (2 * N - 1))
+
+    N = len(t_cc)
+    dnu = 1.0 / (N * dt)
+
+    nu = np.linspace(0, (N - 1) * dnu, N)
+    tau = t_cc
+
+    cc_freqs = scipy.fftpack.fftfreq(len(t_cc), d=dt)
+    freqs = scipy.fftpack.fftfreq(len(t), d=dt)
+
+    # Compute the time frequency representation
+    tfs = np.zeros((len(t), len(t)), dtype="complex128")
+
+    threshold = np.abs(s1).max() * threshold
+
+    for k in range(len(t)):
+        # Window the signals
+        w = utils.gaussian_window(t - tau[k], width)
+        f1 = w * s1
+        f2 = w * s2
+
+        if min(np.abs(f1).max(), np.abs(f2).max()) < threshold:
+            continue
+
+        cc = utils.cross_correlation(f2, f1)
+        tfs[k, :] = scipy.interpolate.interp1d(
+            cc_freqs, scipy.fftpack.fft(cc)
+        )(freqs)
+    tfs *= dt / np.sqrt(2.0 * np.pi)
+
+    return tau, nu, tfs
+
+
+def itfa(tau, tfs, width, threshold=1e-2):
+    N = len(tau)
+    dt = tau[1] - tau[0]
+
+    threshold = np.abs(tfs).max() * threshold
+
+    # inverse fft
+    I = np.zeros((N, N), dtype="complex128")
+
+    # IFFT and scaling.
+    for k in range(N):
+        if np.abs(tfs[k, :]).max() < threshold:
+            continue
+        I[k, :] = scipy.fftpack.ifft(tfs[k, :])
+    I *= 2.0 * np.pi / dt
+
+    # time integration
+    s = np.zeros(N, dtype="complex128")
+
+    for k in range(N):
+        f = utils.gaussian_window(tau[k] - tau, width) * I[:, k].transpose()
+        s[k] = np.sum(f) * dt
+    s *= dt / np.sqrt(2.0 * np.pi)
+
+    return s, tau, I
diff --git a/pyadjoint/timefrequency_utils/utils.py b/pyadjoint/timefrequency_utils/utils.py
new file mode 100644
index 0000000..4e39853
--- /dev/null
+++ b/pyadjoint/timefrequency_utils/utils.py
@@ -0,0 +1,242 @@
+import inspect 
+import matplotlib.pyplot as plt
+import os
+import numpy as np
+
+import obspy
+
+
+EXAMPLE_DATA_PDIFF = (800, 900)
+EXAMPLE_DATA_SDIFF = (1500, 1600)
+
+
+def window_trace(trace, window, taper, taper_ratio, taper_type):
+    """
+    Helper function to taper a window within a data trace.
+
+    This function modifies the passed trace object in-place.
+
+    :param trace: The trace to be tapered.
+    :type trace: :class:`obspy.core.trace.Trace`
+    :param window: Tuples with UCTDateTime objects for start and end time
+        and potentially a weight as well
+    :type window: Tuple with UCTDateTime objects
+    :param taper: True if you want to apply tapering
+    :type taper: binary
+    :param taper_percentage: Decimal percentage of taper at one end (ranging
+        from ``0.0`` (0%) to ``0.5`` (50%)).
+    :type taper_percentage: float
+    :param taper_type: The taper type, supports anything
+        :meth:`obspy.core.trace.Trace.taper` can use.
+    :type taper_type: str
+
+    Any additional keyword arguments are passed to the
+    :meth:`obspy.core.trace.Trace.taper` method.
+    """
+    s, e = trace.stats.starttime, trace.stats.endtime
+   # print(window)
+   # print(trace)
+    # print(s+window[0], s+window[1])
+    # print(taper)
+    trace.trim(s+window[0], s+window[1])
+    # print(trace)
+   # print(len(trace))
+    # print(taper_ratio,taper_type)
+    if taper:
+        trace.taper(max_percentage=taper_ratio, type=taper_type, side='both')
+   # print(len(trace))
+    #
+    trace.trim(s, e, pad=True, fill_value=0.0)
+    # Enable method chaining.
+    return trace
+
+
+def get_example_data():
+    """
+    Helper function returning example data for SalvusMisft.
+
+    :returns: Tuple of observed and synthetic streams
+    :rtype: tuple of :class:`obspy.core.stream.Stream` objects
+
+    .. rubric:: Example
+
+    >>> from lasif.tools.adjoint.utils import get_example_data
+    >>> observed, synthetic = get_example_data()
+    >>> print(observed)  # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
+    3 Trace(s) in Stream:
+    SY.DBO.S3.MXR | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    SY.DBO.S3.MXT | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    SY.DBO.S3.MXZ | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    >>> print(synthetic)  # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
+    3 Trace(s) in Stream:
+    SY.DBO..LXR   | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    SY.DBO..LXT   | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    SY.DBO..LXZ   | 2014-11-15T02:31:50.259999Z - ... | 1.0 Hz, 3600 samples
+    """
+    path = os.path.join(
+        os.path.dirname(inspect.getfile(inspect.currentframe())),
+        "example_data",
+    )
+    observed = obspy.read(os.path.join(path, "observed_processed.mseed"))
+    observed.sort()
+    synthetic = obspy.read(os.path.join(path, "synthetic_processed.mseed"))
+    synthetic.sort()
+
+    return observed, synthetic
+
+
+def generic_adjoint_source_plot(
+    observed, synthetic, time, adjoint_source, misfit, adjoint_source_name
+):
+    """
+    Generic plotting function for adjoint sources and data.
+
+    Many types of adjoint sources can be represented in the same manner.
+    This is a convenience function that can be called by different
+    the implementations for different adjoint sources.
+
+    :param observed: The observed data, windowed.
+    :type observed: numpy array
+    :param synthetic: The synthetic data, windowed.
+    :type synthetic: numpy array
+    :param stats: Stats from obspy trace
+    :type stats: Dictionary
+    :param adjoint_source: The adjoint source.
+    :type adjoint_source: `numpy.ndarray`
+    :param misfit: The associated misfit value.
+    :float misfit: misfit value
+    :param adjoint_source_name: The name of the adjoint source.
+    :type adjoint_source_name: str
+    """
+
+    plt.subplot(211)
+    plt.plot(time, observed, color="0.2", label="Observed", lw=2)
+    plt.plot(time, synthetic, color="#bb474f", label="Synthetic", lw=2)
+    plt.grid()
+    plt.legend(fancybox=True, framealpha=0.5)
+
+    plt.subplot(212)
+    plt.plot(
+        time, adjoint_source, color="#2f8d5b", lw=2, label="Adjoint Source"
+    )
+    plt.grid()
+    plt.legend(fancybox=True, framealpha=0.5)
+
+    plt.suptitle(
+        "%s Adjoint Source with a Misfit of %.3g"
+        % (adjoint_source_name, misfit)
+    )
+
+
+def matlab_range(start, stop, step):
+    """
+    Simple function emulating the behaviour of Matlab's colon notation.
+
+    This is very similar to np.arange(), except that the endpoint is included
+    if it would be the logical next sample. Useful for translating Matlab code
+    to Python.
+    """
+    # Some tolerance
+    if (abs(stop - start) / step) % 1 < 1e-7:
+        return np.linspace(
+            start, stop, int(round((stop - start) / step)) + 1, endpoint=True
+        )
+    return np.arange(start, stop, step)
+
+
+def get_dispersed_wavetrain(
+    dw=0.001,
+    distance=1500.0,
+    t_min=0,
+    t_max=900,
+    a=4,
+    b=1,
+    c=1,
+    body_wave_factor=0.01,
+    body_wave_freq_scale=0.5,
+    dt=1.0,
+):
+    """
+    :type dw: float, optional
+    :param dw: Angular frequency spacing. Defaults to 1E-3.
+    :type distance: float, optional
+    :param distance: The event-receiver distance in kilometer. Defaults to
+        1500.
+    :type t_min: float, optional
+    :param t_min: The start time of the returned trace relative to the event
+        origin in seconds. Defaults to 0.
+    :type t_max: float, optional
+    :param t_max: The end time of the returned trace relative to the event
+        origin in seconds. Defaults to 900.
+    :type a: float, optional
+    :param a: Offset of dispersion curve. Defaults to 4.
+    :type b: float, optional
+    :param b: Linear factor of the dispersion curve. Defaults to 1.
+    :type c: float, optional
+    :param c: Quadratic factor of the dispersion curve. Defaults to 1.
+    :type body_wave_factor: float, optional
+    :param body_wave_factor: The factor of the body waves. Defaults to 0.01.
+    :type body_wave_freq_scale: float, optional
+    :param body_wave_freq_scale:  Determines the frequency of the body waves.
+        Defaults to 0.5
+    :returns: The time array t and the displacement array u.
+    :rtype: Tuple of two numpy arrays
+    """
+    # Time and frequency axes
+    w_min = 2.0 * np.pi / 50.0
+    w_max = 2.0 * np.pi / 10.0
+    w = matlab_range(w_min, w_max, dw)
+    t = matlab_range(t_min, t_max, dt)
+
+    # Define the dispersion curves.
+    c = a - b * w - c * w ** 2
+
+    # Time integration
+    u = np.zeros(len(t))
+
+    for _i in range(len(t)):
+        u[_i] = np.sum(w * np.cos(w * t[_i] - w * distance / c) * dw)
+
+    # Add body waves
+    u += (
+        body_wave_factor
+        * np.sin(body_wave_freq_scale * t)
+        * np.exp(-((t - 250) ** 2) / 500.0)
+    )
+
+    return t, u
+
+
+def cross_correlation(f, g):
+    """
+    Computes a cross correlation similar to numpy's "full" correlation, except
+    shifted indices.
+
+    :type f: numpy array
+    :param f: function 1
+    :type g: numpy array
+    :param g: function 1
+    """
+    cc = np.correlate(f, g, mode="full")
+    N = len(cc)
+    cc_new = np.zeros(N)
+
+    cc_new[0 : (N + 1) // 2] = cc[(N + 1) // 2 - 1 : N]
+    cc_new[(N + 1) // 2 : N] = cc[0 : (N + 1) // 2 - 1]
+    return cc_new
+
+
+def gaussian_window(y, width):
+    """
+    Returns a simple gaussian window along a given axis.
+
+    :type y: numpy array
+    :param y: The values at which to compute the window.
+    :param width: float
+    :param width: variance = (width ^ 2) / 2
+    """
+    return (
+        1.0
+        / (np.pi * width ** 2) ** (0.25)
+        * np.exp(-0.5 * y ** 2 / width ** 2)
+    )
diff --git a/pyadjoint/utils/signal.py b/pyadjoint/utils/signal.py
index 1a90791..0afea80 100644
--- a/pyadjoint/utils/signal.py
+++ b/pyadjoint/utils/signal.py
@@ -270,5 +270,57 @@ def process_cycle_skipping(phi_w, nfreq_max, nfreq_min, wvec, phase_step=1.5):
 
     return phi_w
 
+def matlab_range(start, stop, step):
+    """
+    Simple function emulating the behaviour of Matlab's colon notation.
+
+    This is very similar to np.arange(), except that the endpoint is included
+    if it would be the logical next sample. Useful for translating Matlab code
+    to Python.
+    """
+    # Some tolerance
+    if (abs(stop - start) / step) % 1 < 1e-7:
+        return np.linspace(
+            start, stop, int(round((stop - start) / step)) + 1, endpoint=True
+        )
+    return np.arange(start, stop, step)
+
+def window_trace(trace, window, taper, taper_ratio, taper_type):
+    """
+    Helper function to taper a window within a data trace.
 
+    This function modifies the passed trace object in-place.
+
+    :param trace: The trace to be tapered.
+    :type trace: :class:`obspy.core.trace.Trace`
+    :param window: Tuples with UCTDateTime objects for start and end time
+        and potentially a weight as well
+    :type window: Tuple with UCTDateTime objects
+    :param taper: True if you want to apply tapering
+    :type taper: binary
+    :param taper_percentage: Decimal percentage of taper at one end (ranging
+        from ``0.0`` (0%) to ``0.5`` (50%)).
+    :type taper_percentage: float
+    :param taper_type: The taper type, supports anything
+        :meth:`obspy.core.trace.Trace.taper` can use.
+    :type taper_type: str
 
+    Any additional keyword arguments are passed to the
+    :meth:`obspy.core.trace.Trace.taper` method.
+    """
+    s, e = trace.stats.starttime, trace.stats.endtime
+   # print(window)
+   # print(trace)
+    # print(s+window[0], s+window[1])
+    # print(taper)
+    trace.trim(s + window[0], s + window[1])
+    # print(trace)
+    # print(len(trace))
+    # print(taper_ratio,taper_type)
+    if taper:
+        trace.taper(max_percentage=taper_ratio, type=taper_type, side='both')
+   # print(len(trace))
+    #
+    trace.trim(s, e, pad=True, fill_value=0.0)
+    # Enable method chaining.
+    return trace