function [f_est J] = EMexample1 (f, N, sigma2)

% test with 0.038 0.082
% test with 0.036 0.084

% CREATE x
p = length(f);
x = sqrt(sigma2)*randn(1,N);
for i = 1:p
    x = x + cos(2*pi*f(i)*[0:N-1]);
end;

% f values to try for grid search
f_int = 0:0.0001:0.5;

% initial guess - for 2 cosine example
f_est = [0.04 0.08];
% % randomize initial guess
% f_est = f_est + (rand(1,2)-0.5)*0.06;

f_est

cnt = 0;
change = 1;
while (change>0 && cnt < 20 )
    % E-Step
    temp = zeros(p,N);
    for i = 1:p
        temp(i,:) = cos(2*pi*f_est(i)*[0:N-1]);
    end;
    for i = 1:p
        y(i,:) = temp(i,:) + (1/p)*(x - sum(temp,1));
    end;
    figure(1);clf;
    subplot(3,1,1);hold on;
    plot(x,'k-');
    plot(sum(temp,1),'r-');
    subplot(3,1,2);plot(temp');
    subplot(3,1,3);plot(x - sum(temp,1));
    pause(0.1);
    
    % M-step
    J = 0;
    for i = 1:p
        temp = zeros(1,length(f_int));
        for k = 1:length(f_int)
            temp(k) = sum(y(i,:).*cos(2*pi*f_int(k)*[0:N-1]));
        end;
        [maxJ maxIndex] = max(temp);
        J = J + maxJ;
        f_new(i) = f_int(maxIndex);
    end;
    change = sum(abs(f_new-f_est))/p;
    f_est = f_new
    cnt = cnt + 1;
end;
J