hann_func = @(n,N) 0.5*(1 - cos(2*pi*n/(N-1)));
hann_d_func = @(n,N,fs) (pi*fs/(N-1))*sin(2*pi*n/(N-1));
fs = 44100; N = 8192; f1 = 150; f2 = 420;
t_orig = 0:(1/fs):2.5;
x = cos(2*pi*f1*t_orig);
win = hann_func(1:N,N);
win_d = hann_d_func(1:N, N, fs);
[X, f, t] = stft(x,fs,'Window',win,'OverlapLength',N/2,'fftLength',N);
[Xd, ~, ~] = stft(x,fs,'Window',win_d,'OverlapLength',N/2,'fftLength',N);
X_conj_over_norm = conj(X)./ (sum(abs(X)).^2);
Xdphase_dt = -imag(Xd .* X_conj_over_norm);
reass_f_X = f + Xdphase_dt/(2*pi);