function [dfeat] = diag_feat_T_F(feat, T, F)

% Calculates delta features along time and frequency.
% Special considerations are done for values along edges

[Nl, Ml] = size(feat);

gap = T * ones(1, Ml);
gap(1:T) = [0:(T-1)];
gap((Ml-T+1):Ml) = [(T-1):-1:0];
Ldelta_ids = [1:Ml] +gap;
Rdelta_ids = [1:Ml] -gap;

gap = F * ones(Nl,1);
gap(1:F) = [0:(F-1)];
gap((Nl-F+1):Nl) = [(F-1):-1:0];

Ldiag_ids = [1:Nl]' -gap;
Rdiag_ids = [1:Nl]' +gap;

dfeat = feat(Ldiag_ids, Ldelta_ids) - feat(Rdiag_ids, Rdelta_ids);