% Note: n must be odd.
B = zeros(n,n,n,2);
c = zeros(n,1);
r = zeros(n,1);
q = ceil(n/2);
for k=1:q
    c(1:q) = P(q:n,q+k-1);
    r(1:q) = P(q:-1:1,q+k-1);
    B(:,:,k,1) = toeplitz(c,r);
    c(1:q) = P(q:n,q+1-k);
    r(1:q) = P(q:-1:1,q+1-k);
    B(:,:,k,2) = toeplitz(c,r);
end

A = zeros(n^2,n^2);
for i=1:n
    A( (i-1)*n+1:i*n , (i-1)*n+1:i*n ) = B(:,:,1,1);
    for j=1:i-1
        A( (i-1)*n+1:i*n , (j-1)*n+1:j*n ) = B(:,:,i-j+1,1);
        A( (j-1)*n+1:j*n , (i-1)*n+1:i*n ) = B(:,:,i-j+1,2);
    end
end