covariance of weighted multidimensional samples

I have a multidimensional weighted sample set, e.g.

 D=2; % dimension
 N=100; % number of samples
 x=randn(D,N); % samples
 w=ones(1,N)/N; % corresponding weights

I would like to find the weighted covariance of this sample set. To obtain this, I first computed the weighted mean using the formula \mu = \sum_{i=1}^{N} w_{i} x_{i} as

 mu=sum(bsxfun(@times,w,x),2);

Then I need to find the covariance according to the formula \Sigma = \sum_{i=1}^{N} w_{i} (x_{i} - \mu)*(x_{i} - \mu)'. My current code is this:

    ct=bsxfun(@minus,particles,mu);
    P=zeros(Dx,Dx,N);
    for n=1:N,
        P(:,:,n)=w(n)*(ct(:,n)*ct(:,n)');
    end
    Sigma=sum(P,3);

What is the computationally best coding procedure to calculate this covariance? Suggestions appreciated, thanks.