STFT求解瞬时频率matlab代码

以下是MATLAB代码，用于计算信号的STFT及瞬时频率：

% 定义信号
fs = 1000;  % 采样率
t = 0:1/fs:1;  % 时间序列
f1 = 10;  % 信号频率1
f2 = 50;  % 信号频率2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);  % 信号

% STFT计算
winlen = 0.05;  % 窗口长度（秒）
nfft = 1024;  % FFT点数
overlap = 0.5;  % 重叠率
win = hamming(round(winlen*fs));  % 窗口函数
[S, f, t] = spectrogram(x, win, round(overlap*winlen*fs), nfft, fs);
S = abs(S);  % 取幅值

% 瞬时频率计算
df = fs/nfft;  % 频率分辨率
f0 = zeros(size(S,2),size(S,1));  % 初始化瞬时频率矩阵
for i = 1:size(S,2)
    for j = 2:size(S,1)-1
        f0(i,j) = (arg(S(j+1,i))-arg(S(j-1,i)))/(2*pi*df*winlen);
    end
end

% 画图
figure;
subplot(2,1,1);
plot(t, x);
xlabel('Time (s)');
ylabel('Amplitude');
title('Signal');
subplot(2,1,2);
imagesc(t, f, f0);
set(gca,'YDir','normal');
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Instantaneous frequency');

这段代码首先定义了一个包含两个正弦信号的信号，并利用MATLAB的`spectrogram`函数计算了其STFT。然后根据上面提到的公式，计算了瞬时频率，并将结果绘制成图像。