function Phi = BCH_CS_Matrix(n , k)
% Phi = BCH_CS_Matrix(n , k)
%
% This function constructs an RIP fulfiling bipolar matrix which is based
% on the BCH-based design of the following paper:
%
% [*] A. Amini and F. Marvasti, ?Deterministic Construction of Binary,
% Bipolar and Ternary Compressed Sensing Matrices," IEEE Trans. on
% Inform. Theory, vol. 57, no. 4, pp. 2360-2370, April 2011.
%
% The generated matrix is an m*n column-normalized bipolar matrix with a
% coherence value less than 1/(k-1) which is a sufficient condition for
% establishing the RIP of order "k" and also perfect reconstrcution of
% (k/2)-sparse vectors. In fact, "k", "n" and "m" are the standard notations
% in compressed sensing.
%
% The function is written by Arash Amini
% parameters
RootPower = 1; % choosing the primitive root of the field by
% setting a power for the default primitive root (2) in MATLAB
% determining \tilde{m} for the requested n
i_alg = ceil(log2(k));
tildem = i_alg - 1;
hx = [1];
while 2 ^ (length(hx) - 2) < n
tildem = tildem + 1;
[hx , gx , hxtext] = PolyBCH(tildem , i_alg , RootPower);
end
m = 2 ^ tildem - 1;
if 2 ^ (length(hx) - 2) > n
disp(['!!!! Note: The same number of samples can be used (the same "k") for n <= ' , num2str(2 ^ (length(hx) - 2)) , ' !!!!'])
end
% modifing the code for the even parity
gx = mod([gx 0] + [0 gx] , 2);
gx = double(gx);
% generating the matrix
Phi = zeros(m , n);
for aa = 0 : n - 1
% aa
uncoded = de2bi(aa , 'left-msb' , length(hx) - 2);
Phi(: , aa + 1) = mod(conv(uncoded , gx) , 2).';
end
% Bipolar representation and normalization
Phi(Phi == 0) = -1;
Phi = Phi / sqrt(m);
function [hx , gx , hxtext] = PolyBCH(tildem , i , RootPower)
% [hx , gx , hxtext] = PolyBCH(tildem , i , RootPower)
%
% This function finds the binary parity check ("hx") and code generating
% ("gx") polynomials of the BCH code used for CS purposes.
%
% The roots of these binary polynomials belong to GF(2^tildem); moreover,
% the roots of h(x) are restricted to lie in the set
% {1, alpha, alpha^2, ..., alpha^( 2^(tildem-1) + 2^(tildem - i) )}
% where "alpha" is one of the primitive roots of the field.
%
% The input "RootPower" is an optional input which specifies the primitive
% root. The default value for the root is '2'; however, if "RootPower" is
% given, this root is 2^RootPower.
%
% Note: "RootPower" should be prime relative to 2^tildem - 1; when this
% condition is violated, an error mesage will appear.
%
% Remark: The outputs "hx" and "gx" are binary row vectors which show the
% polynomial coefficients with the descending order with respect
% to the powers of 'x'; i.e., the last elements are constant
% terms of the polynomial. In addition, "hxtext" is a string
% which shows the parity check polynomial.
%
% The function is written by Arash Amini
% setting the primitive root
if nargin == 2
alpha = gf(2 , tildem);
else
if gcd(RootPower , 2^tildem - 1) == 1
alpha = gf(2 , tildem) ^ RootPower;
else
disp(['Error: The input ''RootPower'' should be prime to 2^tildem - 1 which is ' , num2str(2^tildem - 1) , ' in this setting'])
return
end
end
% finding the root powers of the parity check polynomial ( h(x) )
Hbin = binHseq(tildem , i);
Hseq = bi2de(Hbin , 'left-msb');
% evaluating the parity check polynomial by its roots
Hpoly = gf(1 , tildem);
for bb = 1 : length(Hseq)
Hpoly = conv(Hpoly , [alpha^0 alpha^(Hseq(bb))]);
end
hx = uint8(Hpoly.x);
deg_h = length(hx) - 1;
% producing the string that represents the parity check polynomial
hxtext = [' h(x) = 1'];
for bb = deg_h : -1 : 1
if hx(bb) == 1
hxtext = [hxtext , ' + x^' , num2str(length(Hseq) - bb + 1)];
end
end
% finding g(x) by dividing x^(2^tildem - 1) + 1 by h(x)
BigPoly = uint8([1 zeros(1 , 2^tildem - 2) 1]);
gx = uint8(zeros(1 , 2^tildem - deg_h));
for bb = 1 : length(gx)
if BigPoly(bb) == 1
BigPoly(bb : bb + deg_h) = mod(BigPoly(bb : bb + deg_h) + hx , 2);
gx(bb) = 1;
end
end
function Hbin = binHseq(tildem , i)
% Hbin = binHseq(tildem , i)
%
% This function finds all binary sequences of length "tildem" such that 1s
% are circularly spaced with at least $i$ zeros. These binary sequences are
% put together as rows of a matrix ("H")
%
% The function is written by Arash Amini
Allseq = de2bi([0 : 2^tildem - 1] , 'left-msb');
isOK = zeros(1 , 2^tildem);
AllInd = 1 : tildem;
for aa = 1 : 2^tildem
seq = Allseq(aa , :);
loc1 = AllInd(seq == 1);
if length(loc1) > 1
diffs = loc1(2 : end) - loc1(1 : end-1); % spacing between the 1s in the middle
diffs = [loc1(1) + tildem - loc1(end) diffs]; % including the circular spacing between
% the first and the last 1
if min(diffs) > i % whether there exists at least "i" zeros in between
isOK(aa) = 1;
end
else
isOK(aa)= 1;
end
end
Hbin = Allseq(isOK == 1 , :);