This is as much a request for comment and critique as anything else. I am looking at tracking how a system perturbs an input that is a (real) sinoid of varying frequency.
%% generate frequency-modulated sine fs = 48e3; N = fs/2;%2*fs; t = (0:1/fs:((N-1)/fs))'; fc = 4e3; fmod = 2; modulation = cos(2*pi*fmod*t); mod_depth = 3.5e3; y_fwd = exp(+1j*2*pi*(fc.*t+mod_depth*cumsum(modulation)/fs)); y_bwd = exp(-1j*2*pi*(fc.*t+mod_depth*cumsum(modulation)/fs)); y = real(y_fwd); %% introduce time-base disturbance (discontinuity) z = y([1:end/2 (end/2+1:end)-10]);
My goto method for aligning two sequences would be the cross-correlation, plotting here the real component of its output (as y_bwd is complex):
[xk, xklags] = xcorr(z, y_bwd);
Zooming in on the central (0-lag) area, we see an immediate peculiarity: rather than a central maximum at the symmetry center (lag=5), there are two off-center peaks (lag = 0 and 10). This is of course because we have introduced a timeline disontinuity in z.
A tempting direction would be:
%% use Hilbert/analytic signal instfrq_z = fs/(2*pi)*diff(unwrap(angle(hilbert(z))));
As a low complexity alternative, we might count the distance between zero crossings. Note that we are here smoothing a function that is sampled 2 times per (local) cycle, rather than at the uniform samplerate.
%% zero crossings w_z = z(2:end).*z(1:end-1); p_z = (find(w_z<0)); freq_from_zc = fs./(2*diff(p_z)); h_zc = fir1(128-1, 1/7); freq_from_zc_sm = conv(freq_from_zc, h_zc, 'same'); p_z_loc = (p_z(1:end-1) + p_z(2:end))/2;
I was wondering what other solutions one might come up with, and how this relates to traditional tools from electronic engineering/radio. For tracking a slowly varying sinoid, the phase-locked loop (PLL) is a common tool, e.g. employing a feedbackloop of an rf mixer (multiplier) between an input signal and a reference, a smoothing filter and a voltage controlled variable frequency oscillator: https://en.wikipedia.org/wiki/Phase-locked_loop
%% frequency mixer w = z.*y_bwd;
When mixing with the conjugate of the chirp itself, one would ideally expect the filtered output to be a constant. A lowpass filter should leave us with the DC component:
%% frequency mixer, sideband filtering h = fir1(1024-1, (fc-mod_depth)*(2/2)/(fs/2)); w_sm = conv(w, h, 'same');
We can see that the two waveforms are "in lock" for the first part (mean(cos^2(0:x)) tends towards +0.5), while for the latter part, there is a slowly modulated phase.
What else is calculated from z(i)*y(j), i=j? Thats right, the sample correlation matrix (equivalent to the covariance matrix with no mean subtraction and no normalization). At this point I had to limit the 'N' variable to 0.5 seconds due to memory constraints.
corrmat = single(z).'.*single(y_bwd);
The diagonal of the corrmat is basically a Frequency mixer. But the sum of each diagonal also (collectively) forms the cross correlation. Rather than working on the diagonals of that matrix, I prefer to mangle it so each diagonal is now a row in a new, padded matrix:
[r, c] = size(corrmat); corrmat_twist = [zeros([r-1 c], 'single'); corrmat]; rp = 2*r-1; diag_row_idx = (mod([(0:rp-1)'+(c-1:-1:0)], rp)+1)+(0:rp:(c-1)*rp); corrmat_twist = corrmat_twist(diag_row_idx); xk = sum(corrmat_twist,2); corrmat_twist_sm = conv2(corrmat_twist, h, 'same');
Like with the frequency mixer, I'll use the same lowpass filter to smooth my recently "twisted" diagonals along time:
corrmat_twist_sm = conv2(corrmat_twist, h, 'same');
The (real) raw "twisted" xcorr matrix to the left, the (real) sum of each row in the middle (being equivalent to regular cross-correlation), and the "smoothed" xcorr matrix to the right, being that it is a sort of an intermediate between the first two (finite length averaging, rather than full row length averaging).
The idea that brought me here in the first place was that as the global cross-correlation is distorted due to the time-base shift of the latter part, partial cross-correlations in the first or the latter half should be basically perfect (as long as the shift is sufficiently small). I.e. "local calculation" of the cross correlation seems to make sense in this case.
Picking the central few rows of the smoothed twisted xcorr matrix (equivalent with a few diagonals above and below the main diagonal), we see that the magnitude (top) contains few clues in this case. The phase angle is interesting though. Fo the first (undistorted) half, the line consistent with the main diagonal (no offset) is a constant. For the latter part, an offset consistent with 10 samples delay is a constant.
Can this be used for anything? Have I arrived at some for of "PLL" that would be forward controlled rather than a feedback loop?