function [Z1,Z2] = circ_cov_sample(c1)
% C.E. Powell, MATH64082
% Computes two samples Z1,Z2 from N(0, C_N) using DFT method if C_N 
% is a CIRCULANT valid covariance matrix.  
% input: 
%        c1  - first column vector of C_N, of length N
% outputs:
%        Z1, Z2 - two samples from N(0,C_N):

N=length(c1); xi=randn(N,2)*[1; sqrt(-1)];
% eigenvalues of C_N;
d=ifft(c1,'symmetric')*N;
% check minimum eigenvalue 
min(d)
% multiplication with Fourier matrix W
Z=fft((d.^0.5).*xi)/sqrt(N);
% Take real and imag parts;
Z1=real(Z); Z2=imag(Z);
 
