function [A] = generate_matrix(M, N, ensemble, p)
% =========================================================================
%                       Measerement matrix generator
% =========================================================================
% INPUT ARGUMENTS:
% M                         Number of measurements (number of rows).
% N                         Size of sparse vector (number of columns).
% ensemble                  Ensemble type of measurement matrix. Possible
%                           values are:
%                           -'Gaussian': creates a MxN measurement matrix with 
%                           elements drawn from normal distribution N(0,1).%                           
%                           -'Bernoulli': creates a MxN measurement matrix with 
%                           elements drawn from Bernoulli distribution (1/2,1/2).
%                           -'pBernoulli': creates a MxN measurement matrix with 
%                           elements drawn from Bernoulli distribution (p,1-p).
%                           Parameter of Bernoulli distribution.
%                           -'sparseGaussian': creates a MxN sparse measurement 
%                           matrix with elements drawn from normal distribution N(0,1).
% =========================================================================
% OUTPUT ARGUMENTS:
% A                     MxN measurement matrix with normalized columns.
% =========================================================================
% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
% =========================================================================
if nargin < 3
    ensemble = 'Gaussian';
end;

if nargin < 4
    p = 0.5;
end;

switch ensemble
    case 'Gaussian'
        A = randn(M,N);         % Standard normal distribution                
        for i = 1:N             % Normalize columns
            A(:,i) = A(:,i)/norm(A(:,i));
        end;
    case 'Bernoulli'
        A = (-1).^round(rand(M,N));     % Bernoulli ~ (1/2, 1/2) distribution
        for i = 1:N             % Normalize columns
            A(:,i) = A(:,i)/norm(A(:,i));
        end;
    case 'pBernoulli'
        A = (-1).^(rand(M,N) > p);   % Bernoulli ~ (p, 1-p) distribution     
        for i = 1:N             % Normalize columns
            A(:,i) = A(:,i)/norm(A(:,i));
        end;
    case 'sparseGaussian'
        leftd = 8;
        A = zeros(M,N);
        for i = 1:N
            ind = randperm(M);
            A(ind(1:leftd),i)=1/leftd;
        end
        for i = 1:N             % Normalize columns
            A(:,i) = A(:,i)/norm(A(:,i));
        end;
end;