Определение зависящей от времени частоты с использованием скользящего окна FFT

Question

У меня есть инструмент, который выдает примерно синусоидальные данные, но с частотой, слегка изменяющейся во времени. Я использую MATLAB для создания прототипа некоторого кода для характеристики зависимости от времени, но я сталкиваюсь с некоторыми проблемами.

Я генерирую идеализированное приближение моих данных, I(t) = sin (2 pi f(t) t), с переменной f(t), но в настоящее время проверяется как линейная или квадратичная. Затем я реализую скользящее окно Хэмминга (шириной w), чтобы сгенерировать набор преобразований Фурье F[I(t), t'], соответствующих точкам данных в I(t) и каждому F[I(t), t'] соответствует гауссову для более точного определения местоположения пика.

Мой текущий код MATLAB:

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

Мой мыслительный процесс состоит в том, что местоположение пика в каждом F[I(t), t'] должно быть близким приближением частоты в центре окна, которое использовалось для его создания. Однако, похоже, что это не так, экспериментально.

Я имел некоторый успех, используя дискретный анализ Фурье для инженерных задач в прошлом, но я только делал курсовую работу по непрерывным преобразованиям Фурье - так что может быть что-то очевидное, что я упускаю. Кроме того, это мой первый вопрос о StackExchange, так что конструктивная критика приветствуется.

Answer 1

Вот и получается, что моей проблемой было плохое понимание математики функции синуса. Я предположил, что частота волны была равна тому, что было умножено на переменную времени (например, f в грехе (футах)). Однако оказывается, что частота фактически определяется производной всего аргумента функции синуса - скоростью изменения фазы.

Для константы f оба определения равны, поскольку d (ft)/dt = f. Но, скажем, f (t) = sin (t):

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

Частота изменяется как функция, сильно отличающаяся от f (t). Изменение определения функции на следующее решило мою проблему:

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