使用滑动窗口 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

所以事实证明我的问题是对正弦函数的数学理解不足。我假设波的频率等于乘以时间变量的值（例如 sin(ft) 中的 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