function DTFTgram_cs_spectrum = DTFT_circ_precon(omegas,Nd)
% 
% Copyright 2013, Daniel Weller, University of Michigan. All rights
% reserved. Email wellerds@umich.edu.

% construct DTFT gram matrix (hermitian p.s.d., Toeplitz)
% tinds = cell(1,length(Nd));
% for id = 1:length(Nd), tinds{id} = -Nd(id)+1:Nd(id)-1; end
tinds = arrayfun(@(n) -n+1:n-1,Nd,'UniformOutput',false);
[tinds{:}] = ndgrid(tinds{:});
tinds = cellfun(@(x) col(x).', tinds, 'UniformOutput', false);
tinds = num2cell(vertcat(tinds{:}),1);

DTFTgram_ts = cellfun(@(x) sum(exp(complex(0,omegas*x))),tinds);

DTFTgram_ts = DTFTgram_ts./prod(Nd); % "unitary" NUFFT
DTFTgram_ts = reshape(DTFTgram_ts,[2*Nd-1,1]);

% compute Frobenius-norm-optimal circulant matrix of same size
DTFTgram_cs = DTFTgram_ts;
for id = 1:length(Nd)
    weight_cs_pos = col(Nd(id):-1:1)./Nd(id);
    weight_cs_neg = col(0:Nd(id)-1)./Nd(id);
    DTFTgram_cs = repmat(weight_cs_pos,[1,2*Nd(id+1:end)-1,Nd(1:id-1)]).*DTFTgram_cs(Nd(id):end,:) ...
        + repmat(weight_cs_neg,[1,2*Nd(id+1:end)-1,Nd(1:id-1)]).*DTFTgram_cs([Nd(id),1:Nd(id)-1],:);
    DTFTgram_cs = shiftdim(DTFTgram_cs,1);
end

% note that circulant matrix is symm. pos. def. whenever DHD_alphaI is, which
% should be always
DTFTgram_cs_spectrum = fftn(DTFTgram_cs);
DTFTgram_cs_spectrum = reale(DTFTgram_cs_spectrum,'warn','imaginary eigenvalue, check symmetry');
DTFTgram_cs_cond = max(DTFTgram_cs_spectrum(:))/min(DTFTgram_cs_spectrum(:));
if DTFTgram_cs_cond < 0
    warn('negative eigenvalue, check positive semi-definiteness');
elseif ~isfinite(DTFTgram_cs_cond)
    fprintf(1,'DTFT gram matrix: zero eigenvalue, requires regularization\n');
end

end
