Determining time-dependent frequency using a sliding-window FFT

Question

I have an instrument which produces roughly sinusoidal data, but with frequency varying slightly in time. I am using MATLAB to prototype some code to characterize the time dependence, but I'm running into some issues.

I am generating an idealized approximation of my data, I(t) = sin(2 pi f(t) t), with f(t) variable but currently tested as linear or quadratic. I then implement a sliding Hamming window (of width w) to generate a set of Fourier transforms F[I(t), t'] corresponding to the data points in I(t), and each F[I(t), t'] is fit with a Gaussian to more precisely determine the peak location.

My current MATLAB code is:

fs = 1000; %Sample frequency (Hz)
tlim = [0,1];
t = (tlim(1)/fs:1/fs:tlim(2)-1/fs)'; %Sample domain (t)
N = numel(t);

f = @(t) 100-30*(t-0.5).^2; %Frequency function (Hz)
I = sin(2*pi*f(t).*t); %Sample function

w = 201; %window width
ww=floor(w/2); %window half-width

for i=0:2:N-w

    %Take the FFT of a portion of I, convolved with a Hamming window
    II = 1/(fs*N)*abs(fft(I((1:w)+i).*hamming(w))).^2;
    II = II(1:floor(numel(II)/2));
    p = (0:fs/w:(fs/2-fs/w))';

    %Find approximate FFT maximum
    [~,maxIx] = max(II);
    maxLoc = p(maxIx);

    %Fit the resulting FFT with a Gaussian function
    gauss = @(c,x) c(1)*exp(-(x-c(2)).^2/(2*c(3)^2));
    op = optimset('Display','off');
    mdl = lsqcurvefit(gauss,[max(II),maxLoc,10],p,II,[],[],op);    

    %Generate diagnostic plots
    subplot(3,1,1);plot(p,II,p,gauss(mdl,p))
    line(f(t(i+ww))*[1,1],ylim,'color','r');

    subplot(3,1,2);plot(t,I);
    line(t(1+i)*[1,1],ylim,'color','r');line(t(w+i)*[1,1],ylim,'color','r')

    subplot(3,1,3);plot(t(i+ww),f(t(i+ww)),'b.',t(i+ww),mdl(2),'r.');
    hold on
    xlim([0,max(t)])
    drawnow
end
hold off

My thought process is that the peak location in each F[I(t), t'] should be a close approximation of the frequency at the center of the window which was used to produce it. However, this does not seem to be the case, experimentally.

I have had some success using discrete Fourier analysis for engineering problems in the past, but I've only done coursework on continuous Fourier transforms--so there may be something obvious that I'm missing. Also, this is my first question on StackExchange, so constructive criticism is welcome.

Answer 1

So it turns out that my problem was a poor understanding of the mathematics of the sine function. I had assumed that the frequency of the wave was equal to whatever was multiplied by the time variable (e.g. the f in sin(ft)). However, it turns out that the frequency is actually defined by the derivative of the entire argument of the sine function--the rate of change of the phase.

For constant f the two definitions are equal, since d(ft)/dt = f. But for, say, f(t) = sin(t):

d(f(t)t)/dt = d(sin(t) t)/dt = t cos(t) + sin(t)

The frequency varies as a function very different from f(t). Changing the function definition to the following fixed my problem:

f = @(t) 100-30*(t-0.5).^2; %Frequency function (Hz)
G = cumsum(f(t))/fs; %Phase function (Hz)
I = sin(2*pi*G); %Sampling function