function [P1,P2,P3]=testme(K,L,N,sigma,numobs)

  P = zeros(L,L);
  step = floor(L/3);
  for (k=1:step)
    P(k,k) = 0.5;
    P(k+step,k+step) = 1;
    P(k+2*step,k+2*step) = 1.5;
  end
  
  c1 = N/K;
  c2 = L/N;
  c3 = L/K;
  
    obsnr=numobs;
    samplecovsum1 = 0;
    samplecovsum2 = 0;
    samplecovsum3 = 0;
    for (thisobs=1:obsnr)
      V = zeros(N,L);
      for (k=1:L) 
        V(:,k) = ((exp(2*pi*j*rand(1))).^(0:(N-1)))';
      end
      V = (1/sqrt(N))*V;
      alpha = randn(L,K);
      nmatrix = sigma*randn(N,K);
      Y = V*sqrt(P)*alpha + nmatrix;
      W = (1/K)*Y*Y';%sample covariance matrix sample
      samplecovsum1 = samplecovsum1 + (1/N)*trace(W);
      samplecovsum2 = samplecovsum2 + (1/N)*trace(W^2);
      samplecovsum3 = samplecovsum3 + (1/N)*trace(W^3);
    end
    W1 = samplecovsum1/obsnr;
    W2 = samplecovsum2/obsnr;
    W3 = samplecovsum3/obsnr;
 
    P1 = (W1 -sigma^2)/c2;
    P2 = ( W2     -     c2*(c2+c3)*P1^2     -     2*(sigma^2)*(c2+c3)*P1     -     (sigma^4)*(1+c1) )/( c2*(1-1/N) );
    P3 = ( W3     -     (1-1/N)*(3*(c2^2)*(1+1/(K^2))+3*c2*c3)*P1*P2  ...
                  -     ( (c2^3)*(1+1/(K^2)) + 3*(c2^2)*c3 + c2*(c3^2) )*(P1^3) ...    
                  -     3*(sigma^2)*( (1+c1)*c2 + c1*(c2^2)/(K*L) )*(1-1/N)*P2 ...     
                  -     3*(sigma^2)*( c1*(c2^3)/(K*L) + c2^2 + c3^2 + 3*c2*c3 )*(P1^2) ...    
                  -     3*(sigma^4)*( c1^2+3*c1+1+1/(K^2))*c2*P1 ...   
                  -     (sigma^6)*(c1^2+3*c1+1+1/K^2) ...
         )/( c2*(1+1/(K^2))*(1-3/N+2/(N^2)) );