Identification of time-dependent signals from a spectrogram

Question

I am trying to trace several signals through time from a spectrogram. The issue is that as time evolves, the signals shift in frequency, so can be difficult to isolate them. (See images below) Is there a python function out there that is capable of tracking signals through time? Note that I do not have access to the raw data prior to any FTs - only a 2D array of the spectrogram itself.

What I have tried (to moderate success:)

1. Slice up the 2D spectrogram by time, and identify peaks in the frequency at each time.
2. Try to follow peak locations across time to capture the evolution of each signal as a function of time.

This approach is patchwork code, but the results look promising, as seen here:

The original spectrogram:

Frequency tracking across time:

What I am looking for is a more robust method to accomplish this. I'm hoping to not have to write my own script but to use some developed tool that is able to preform this data analysis. Does such a python tool exist?

Thanks.

Here is the code I used to make these plots. Fully open to complete overhauls as this is obviously not robust.

import numpy as np
from matplotlib.pyplot import *
from matplotlib.colors import LogNorm

data = np.load('spectrogram.npz')

## ---- trim data ---- ###

it1=1000
it2=1700
if1=0
if2=1900
tvec = data['tvec'][it1:it2]
fvec = data['fvec'][if1:if2]
spect = data['spect'][0][if1:if2,it1:it2]
v = np.linspace(min(spect.flatten()),0.1*max(spect.flatten()),100)

## ---- spectrogram ---- ###

figure()
contourf(tvec,fvec,spect,v,cmap=plt.cm.jet,norm=LogNorm())
xlabel('Time (ms)')
ylabel('Frequency (Hz)')

## ---- find modes ---- ###

# initialize to NaNs
kHz12 = np.full(len(tvec),np.nan)
kHz25 = np.full(len(tvec),np.nan)
kHz50 = np.full(len(tvec),np.nan)
kHz50_f0 = 50E3
kHz70 = np.full(len(tvec),np.nan)
for i in range(it2-it1):
    x = spect[:,i]
    peaks, _ = scipy.signal.find_peaks(x,height=100,prominence=50)
    #scatter([tvec[i]]*len(peaks),fvec[peaks],1,marker='o')

    ### --- scroll through time and keep track of each frequency --- ###
    kHz12_subsel = [y for j,y in zip(peaks,x[peaks]) if fvec[j] > 11E3 and fvec[j] < 15E3]
    if kHz12_subsel: kHz12[i] = fvec[peaks][np.where(x[peaks] == np.max(kHz12_subsel))]

    kHz25[i] = fvec[peaks][np.argmax(x[peaks])]

    kHz50_subsel = [y for j,y in zip(peaks,x[peaks]) if fvec[j] > 45E3 and fvec[j] < 60E3]
    if kHz50_subsel: 
        kHz50[i] = fvec[peaks][np.argmin([abs(f-kHz50_f0) for f in fvec[peaks]])]
        kHz50_f0 = kHz50[i]

    kHz70_subsel = [j for j in peaks if fvec[j] > 60E3 and fvec[j] < 80E3]
    if kHz70_subsel:
        kHz70[i] = np.average(fvec[kHz70_subsel])

## ---- plot modes ---- ###

plot(tvec,kHz12,lw=2,c='k')
plot(tvec,kHz25,lw=2,c='k')
plot(tvec,kHz50,lw=2,c='k')
plot(tvec,kHz70,lw=2,c='k')