我有一个仪器可以产生大致正弦数据,但频率随时间略有变化。我正在使用 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 上的第一个问题,因此欢迎提出建设性的批评。