Mixture of Gaussian code is running very very slow.

Here I am pasting my code. Plz tell me if there is any correction to work it fast.......

% This m-file implements the mixture of Gaussians algorithm for background 
% subtraction.  It may be used free of charge for any purpose (commercial
% or otherwise), as long as the author (Seth Benton) is acknowledged.

clear all

% source = aviread('C:\Video\Source\traffic\san_fran_traffic_30sec_QVGA');
%source = aviread('..\test_video\san_fran_traffic_30sec_QVGA_Cinepak');
[FileName, PathName] = uigetfile('*.avi', 'Pick an avi-file');
source= aviread(strcat(PathName, FileName));

% -----------------------  frame size variables -----------------------

fr = source(1).cdata;           % read in 1st frame as background frame
fr_bw = rgb2gray(fr);     % convert background to greyscale
fr_size = size(fr);             
width = fr_size(2);
height = fr_size(1);
fg = zeros(height, width);
bg_bw = zeros(height, width);

% --------------------- MOG variables -----------------------------------

C = 3;                                  % number of gaussian components (typically 3-5)
M = 3;                                  % number of background components
D = 2.5;                                % positive deviation threshold
alpha = 0.10;                           % learning rate (between 0 and 1) (from paper 0.01)
thresh = 0.25;                          % foreground threshold (0.25 or 0.75 in paper)
sd_init = 6;                            % initial standard deviation (for new components) var = 36 in paper
w = zeros(height,width,C);              % initialize weights array
mean = zeros(height,width,C);           % pixel means
sd = zeros(height,width,C);             % pixel standard deviations
u_diff = zeros(height,width,C);         % difference of each pixel from mean
p = alpha/(1/C);                        % initial p variable (used to update mean and sd)
rank = zeros(1,C);    % rank of components (w/sd)
%cmap = rand(236, 3);
%Mov1=zeros(length(source));             %preallocating the array for Mov1
%Mov2=zeros(length(source));             %preallocating the array for Mov2

% --------------------- initialize component means and weights -----------

pixel_depth = 8;                        % 8-bit resolution
pixel_range = 2^pixel_depth -1;         % pixel range (# of possible values)

for i=1:height
    for j=1:width
        for k=1:C

              mean(i,j,k) = rand*pixel_range;     % means random (0-255)
              w(i,j,k) = 1/C;                     % weights uniformly dist
              sd(i,j,k) = sd_init;                % initialize to sd_init

          end
      end
  end

%--------------------- process frames -----------------------------------

for n = 1:length(source)

      fr = source(n).cdata;       % read in frame
      fr_bw = rgb2gray(fr);       % convert frame to grayscale

      % calculate difference of pixel values from mean
      for m=1:C
          u_diff(:,:,m) = abs(double(fr_bw) - double(mean(:,:,m)));
      end

      % update gaussian components for each pixel
      for i=1:height
          for j=1:width

              match = 0;
              for k=1:C                       
                  if (abs(u_diff(i,j,k)) <= D*sd(i,j,k))       % pixel matches component

                      match = 1;                          % variable to signal component match

                      % update weights, mean, sd, p
                      w(i,j,k) = (1-alpha)*w(i,j,k) + alpha;
                      p = alpha/w(i,j,k);                  
                      mean(i,j,k) = (1-p)*mean(i,j,k) + p*double(fr_bw(i,j));
                      sd(i,j,k) =   sqrt((1-p)*(sd(i,j,k)^2) + p*((double(fr_bw(i,j)) - mean(i,j,k)))^2);
                  else                                    % pixel doesn't match component
                      w(i,j,k) = (1-alpha)*w(i,j,k);      % weight slighly decreases

                  end
              end

              w(i,j,:) = w(i,j,:)./sum(w(i,j,:));

              bg_bw(i,j)=0;
              for k=1:C
                  bg_bw(i,j) = bg_bw(i,j)+ mean(i,j,k)*w(i,j,k);
              end

              % if no components match, create new component
              if (match == 0)
                  [min_w, min_w_index] = min(w(i,j,:));  
                  mean(i,j,min_w_index) = double(fr_bw(i,j));
                  sd(i,j,min_w_index) = sd_init;
              end

              rank = w(i,j,:)./sd(i,j,:);             % calculate component rank
              rank_ind = [1:1:C];

              % sort rank values
              for k=2:C               
                  for m=1:(k-1)

                      if (rank(:,:,k) > rank(:,:,m))                     
                          % swap max values
                          rank_temp = rank(:,:,m);  
                          rank(:,:,m) = rank(:,:,k);
                          rank(:,:,k) = rank_temp;

                          % swap max index values
                          rank_ind_temp = rank_ind(m);  
                          rank_ind(m) = rank_ind(k);
                          rank_ind(k) = rank_ind_temp;    

                      end
                  end
              end

              % calculate foreground
              match = 0;
              k=1;

              fg(i,j) = 0;
              while ((match == 0)&&(k<=m))

                  if (w(i,j,rank_ind(k)) >= thresh)
                      if (abs(u_diff(i,j,rank_ind(k))) <= D*sd(i,j,rank_ind(k)))
                          fg(i,j) = 0;
                          match = 1;
                      else
                          fg(i,j) = fr_bw(i,j);     
                      end
                  end
                  k = k+1;
              end
          end
      end

      figure(1),subplot(3,1,1),imshow(fr)
      subplot(3,1,2),imshow(uint8(bg_bw))
      subplot(3,1,3),imshow(uint8(fg)) 

      Mov1(n)  = im2frame(uint8(fg),gray);           % put frames into movie
      Mov2(n)  = im2frame(uint8(bg_bw),gray);           % put frames into movie

end

movie2avi(Mov1,'mixture_of_gaussians_output','fps',50);           % save movie as avi 
movie2avi(Mov2,'mixture_of_gaussians_background','fps',50);           % save movie as avi

[Walter Roberson]

step #0: get rid of the "clear all".

[Duke Watson]

Yeah i have done that but still the problem of video is same. It is running very slow. You can test this code on your pc

[Ahmed]

@Duke Watson, did you come up with the solution for the speed? please let me know

[Algorithms Analyst]

Apply matlab parallel computing toolbox...

like matlab pool open