I wrote a C function to reconstruct the amlplitude and frequency of a line. It can be added as a block into a real time monitor. The idea is to use it to track in real time the frequency and amplitude of the disk modes, during the ring down. I did some tests and finally managed to get the function to compile and run on the cymac2 (the crackling lab cymac).
The following plot shows a simulation, since I can't run the code on the new cymac and I don't have the disk installed anymore. The top panel gives the amplitude of a decaying line, and the bottom panel the frequency offset from a reference local oscillator (more below). The nominal values are an initial amplitude of 1, frequency of 1109.375 Hz to be compared to a 1109.0 local oscillator. The fitted decay time is 10.005 seconds, to be comapred with 10 seconds nominal. There is some additive gaussian noise, that causes the ring down to be unmesurable after about 70 seconds of data.
This code will be used for real time estimation of the disk modes, once theot frequency has been roughly estimated with FFTs. The estimation of the frequency work remarkably well. In the first 20 seconds the mean value is 0.3747 Hz, with a standard deviation of 1.5 mHz. When the SNR gets worse (between t=30s and t=50s) the mean value is 0.3745 with a standard deviation of 20 mHz.
Because of the way it's built (see below) the code is sensitive to DC offsets, so the input signal must be high-passed.
Details on the implementation
The code is based on demodulation of the input signal with a reference local oscillator thta must have a frequency as close as possible to the line we want to track. The inputs to the block are: the signal to be monitored, sine of the local oscillator, cosine of the local oscillator. The outputs are: amplitude squared of the peak, frequency offset in Hz from the local oscillator.
Here's the math. Let's assume that the signal is
and the local oscillator has a frequency f0:
The code multiply the signal by the two local oscillators and average the result over 65kHz / 8 Hz samples. Therefore we get two output streams at 8 Hz which are
Then the sum of the two squared 8 Hz streams give an estimate of the amplitude squared. The code computes this every second
while the arctangent of the ratio gives a phase that varies linearly with time.
For each of the 8Hz samples the code computes the arctangent (using a home-brewed function based on a lookup table, since we can't import math.h in the RCG). It unwraps it, and then every second fit a line to the unwrapped arctangent, to estimate the frequency offset with respect to the local oscillator.
The C function has some parameters hard coded: the main sampling frequency (65536 Hz), the number of points per second to use for the frequency estimation (8 Hz), the fact that the output is computed every second. The first two parameters can be changed, the third one cannot for the moment being.