I am currently implementing acoustic FSK modulation and demodulation. I am not a signal processing guy…
Since you say you have matched filters, and you mention non-coherent detection, I think you're pretty much of a digital communication person already – the step to being a DSP person is pretty small :)
The fully-fledged synchronizer SDR approach
So, the classical way to do this is to take the outputs of your two filters, (magnitude) square them, low pass filter the result to something above $2 f_\text{symbol,max}$, and during/after filtering decimate (if possible), so that we don't waste operations per second.
I made myself a quick real-value bandpass signal FSK demod (I think if I did this again, I would have replaced the two bandpass filters with one combined frequency translation and one low pass filter, giving me complex baseband, or threw a polyphase filterbank at the problem, but whatever) in GNU Radio with the help of the GNU Radio companion (flow graph file here):
Point is that we probably want to symbol timing synchronize between the divide and the threshold block. We could do afterwards – the classical, "microcontroller & digital hardware design" approach would in fact be Muller & Mueller clock sink and that would work there just as well – but let's not waste the slope information.
So, I'd add a polyphase clock sync here; I must admit I could not formulate this better than Tom, a disciple of harris, did in the documentation of the eponymous GNU Radio block:
Timing synchronizer using polyphase filterbanks.
This block performs timing synchronization for PAM signals by minimizing the derivative of the filtered signal, which in turn maximizes the SNR and minimizes ISI.
This approach works by setting up two filterbanks; one filterbank contains the signal's pulse shaping matched filter (such as a root raised cosine filter), where each branch of the filterbank contains a different phase of the filter. The second filterbank contains the derivatives of the filters in the first filterbank. Thinking of this in the time domain, the first filterbank contains filters that have a sinc shape to them. We want to align the output signal to be sampled at exactly the peak of the sinc shape. The derivative of the sinc contains a zero at the maximum point of the $\text{sinc}$ ($\text{sinc}(0) = 1,\, \text{sinc}'(0) = 0$). Furthermore, the region around the zero point is relatively linear. We make use of this fact to generate the error signal.
If the signal out of the derivative filters is $d_i[n]$ for the $i$th filter, and the output of the matched filter is $x_i[n]$, we calculate the error as:
$$e[n] = \frac{\Re\{x_i[n]\} \cdot \Re\{d_i[n]\} + \Im\{x_i[n]\} \cdot \Im\{d_i[n]\}}{2.0}\text{ .}$$
This equation averages the error in the real and imaginary parts. There are two reasons we multiply by the signal itself. First, if the symbol could be positive or negative going, but we want the error term to always tell us to go in the same direction depending on which side of the zero point we are on. The sign of $x_i[n]$ adjusts the error term to do this. Second, the magnitude of $x_i[n]$ scales the error term depending on the symbol's amplitude, so larger signals give us a stronger error term because we have more confidence in that symbol's value. Using the magnitude of $x_i[n]$ instead of just the sign is especially good for signals with low SNR.
The error signal, $e[n]$, gives us a value proportional to how far away from the zero point we are in the derivative signal. We want to drive this value to zero, so we set up a second order loop. We have two variables for this loop; $d_k$ is the filter number in the filterbank we are on and $d_\text{rate}$ is the rate which we travel through the filters in the steady state. That is, due to the natural clock differences between the transmitter and receiver, $d_\text{rate}$ represents that difference and would traverse the filter phase paths to keep the receiver locked. Thinking of this as a second-order PLL, the $d_\text{rate}$ is the frequency and $d_k$ is the phase. So we update $d_\text{rate}$ and $d_k$ using the standard loop equations based on two error signals, $d_\alpha$ and $d_\beta$. We have these two values set based on each other for a critically damped system, so in the block constructor, we just ask for "gain," which is $d_\alpha$ while $d_\beta$ is equal to $\frac{\text{gain}^2}{4}$.
Realizing this is audio frequency
So, parameterizing this block properly (and probably reducing bpf_decim andlpf_decim to give the synchronizer more headroom to shift phases) you could build a very stable symbol timing recovery, that would be totally overkill for your application 😁
Since you're doing sampling with a sound card on both ends, and since I presume you don't have to deal with Doppler, the symbol rate has a fixed relation to the carrier frequencies (e.g. the transmitting side has a higher carrier with a period of 8 sound card samples, and a symbol always takes let's say 128 samples, so there's a fixed ratio), you can take a short cut:
If you know the ratio between carrier frequency and symbol rate, just estimate one and you get the other for free! In this case, FSK, that's particularly easy: just wait for a period when the e.g. upper band is active, and use a quadrature demodulator (i.e. practically $\frac{d\arctan}{dt}[n]$) to give you an estimate of the carrier frequency after your decimating band pass filter, and use the resulting value as frequency correction factor – you immediately solve the carrier frequency recovery and symbol rate recovery problem at once. I don't even think you'd need timing recovery after that – after a thresholding, if you oversample enough (which, at these rates, isn't a problem), you can just make a majority decision in a window, or shift that window as long as necessary to get a maximum in the window.
Regarding Packet Sync
Well, since you've got a working demodulator that will probably run too much out of sync when there's no signal present, just use a fixed known data sequence to find your packet.
