function fig10_14

%%%%  ICA using s^4 on linear waves: K and the optimal solution
%%%%  random data
tic
% generate data
tmax=10;
NN=75;
theta0=zeros(NN,1);
theta00=zeros(NN,1);
% IR=number of cases
IR=10;
% IRR=the paricular case plotted
IRR=randi(IR);
% m=number of data points used (this is N in text)
m=3;
for ic=1:NN
    %     if m<40
    %         m=m+1;
    %     else
    %         m=m+40;
    %     end
    m=m+2;
    num(ic)=m;
    val=0;
    for ir=1:IR
        t=tmax*rand(m,1);
        w1=pi;
        w2=1.7*pi;
        x=zeros(m,2);

        % mixing coefficients
        M11=2; M12=1;
        M21=3; M22=-2;
        for it=1:m
            S(it,1) = sin(w1*t(it));
            S(it,2) = sin(w2*t(it)-5);
        end

        %%%% center S
        sM=sum(S)/m;
        S=S-sM;

        %%% whiten S
        [UR,SR,VR]=svd(S);
        D=[sqrt(m)/SR(1,1) 0;0 sqrt(m)/SR(2,2)];
        Fac=D*VR';
        for it=1:m
            W1(it)=Fac(1,1)*S(it,1)+Fac(1,2)*S(it,2);
            W2(it)=Fac(2,1)*S(it,1)+Fac(2,2)*S(it,2);
        end

        for it=1:m
            X(it,1)=M11*S(it,1) + M12*S(it,2);
            X(it,2)=M21*S(it,1) + M22*S(it,2);
        end

        xM=sum(X)/m;
        X=X-xM;
        XX=X'*X/m;
        [UX,SX,VX]=svd(XX);
        D5=[1/sqrt(SX(1,1)) 0; 0 1/sqrt(SX(2,2))];
        Y=X*VX*D5;

        ntheta=1000;
        theta=linspace(0,2*pi,ntheta);
        for ia=1:ntheta
            V=[[ cos(theta(ia)) -sin(theta(ia))];[ sin(theta(ia)) cos(theta(ia))]];
            ssk=Y*V';
            kks=ssk.^4/m ;
            suk=sum(kks);
            K(ia)=abs(suk(1)-3)+abs(suk(2)-3);
        end

        iz=0;
        for ia=2:ntheta-1
            if K(ia-1) < K(ia) && K(ia+1) < K(ia)
                iz=iz+1;
                %thetaMM(iz)=theta(ia);
                break
            end
        end

        if ir==IRR
            theta00(ic)=theta(ia);
        end
        val=val+theta(ia);
    end
    theta0(ic)=val/IR;
end

clf
%get(gcf);
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.03,'MB',0.2,'MR',0.01,'ML',0.06,'P',0.02,'S',0.02)

co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

plot(num,theta0,'LineWidth',1.6)
hold on
plot(num,theta00,'--r','LineWidth',2)
axis([0 150 0 pi/2])
xlabel('N-axis')
ylabel('Theta')
set(gca,'ytick',[0 pi/4 pi/2])
yticklabels({'0','\pi/4','\pi/2'})
grid on
box on
set(gca,'FontSize',16,'FontWeight','bold')

toc

%exportgraphics(gcf,'/Users/mark/Desktop/subdata.eps')


















