%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dWf = calcdWf(Wf,ds,dx,hT,K);
global Rx m dsdshT dsdxhT

% calculate the gradient of the objective function
dWf = zeros(ds,dx,hT);
hWf = hermite3(Wf);
for k = 1:K
WfRx = fast_ip(Wf,Rx(:,:,:,k));
Eomega = fast_ip(WfRx,hWf);
Eomega(dsdshT) = 0; % this uses the optimal lambdaS
dWf = dWf + fast_ip(Eomega,WfRx);
end
dWf = 2*dWf;

%% apply the straightforward power normalization
dWf = m.*dWf;

% enforce unit gains of the diagonal filters
dWf(dsdxhT) = 0;