Plot histogram and estimated PDF in Matlab

Key focus: With examples, let’s estimate and plot the probability density function of a random variable using Matlab histogram function.

Generation of random variables with required probability distribution characteristic is of paramount importance in simulating a communication system. Let’s see how we can generate a simple random variable, estimate and plot the probability density function (PDF) from the generated data and then match it with the intended theoretical PDF. Normal random variable is considered here for illustration. Other types of random variables like uniform, Bernoulli, binomial, Chi-squared, Nakagami-m are illustrated in the next section.

Note: If you are inclined towards programming in Python, visit this article

Step 1: Create the random variable

A survey of commonly used fundamental methods to generate a given random variable is given in [1]. For this demonstration, we will consider the normal random variable with the following parameters : – mean and – standard deviation. First generate a vector of randomly distributed random numbers of sufficient length (say 100000) with some valid values for and . There are more than one way to generate this. Some of them are given below.

This article is part of the book
Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

● Method 1: Using the in-built random function (requires statistics toolbox)

mu=0;sigma=1;%mean=0,deviation=1
L=100000; %length of the random vector
R = random('Normal',mu,sigma,L,1);%method 1

●  Method 2: Using randn function that generates normally distributed random numbers having and = 1

mu=0;sigma=1;%mean=0,deviation=1
L=100000; %length of the random vector
R = randn(L,1)*sigma + mu; %method 2

● Method 3: Box-Muller transformation [2] method using rand function that generates uniformly distributed random numbers

 mu=0;sigma=1;%mean=0,deviation=1
L=100000; %length of the random vector
U1 = rand(L,1); %uniformly distributed random numbers U(0,1)
U2 = rand(L,1); %uniformly distributed random numbers U(0,1)
Z = sqrt(-2log(U1)).cos(2piU2);%Standard Normal distribution
R = Z*sigma+mu;%Normal distribution with mean and sigma

Step 2: Plot the estimated histogram

Typically, if we have a vector of random numbers that is drawn from a distribution, we can estimate the PDF using the histogram tool.  Matlab supports two in-built functions to compute and plot histograms:

● hist – introduced before R2006a
● histogram – introduced in R2014b

Which one to use ? Matlab’s help page points that the hist function is not recommended for several reasons and the issue of inconsistency is one among them. The histogram function is the recommended function to use.

Estimate and plot the normalized histogram using the recommended ‘histogram’ function. And for verification, overlay the theoretical PDF for the intended distribution. When using the histogram function to plot the estimated PDF from the generated random data, use ‘pdf’ option for ‘Normalization’ option. Do not use the ‘probability’ option for ‘Normalization’ option, as it will not match the theoretical PDF curve.

histogram(R,'Normalization','pdf'); %plot estimated pdf from the generated data

X = -4:0.1:4; %range of x to compute the theoretical pdf
fx_theory = pdf('Normal',X,mu,sigma); %theoretical normal probability density
hold on; plot(X,fx_theory,'r'); %plot computed theoretical PDF
title('Probability Density Function'); xlabel('values - x'); ylabel('pdf - f(x)'); axis tight;
legend('simulated','theory');
Estimated PDF (using histogram function) and the theoretical PDF
Estimated PDF (using histogram function) and the theoretical PDF

However, if you do not have Matlab version that was released before R2014b, use the ‘hist’ function and get the histogram frequency counts () and the bin-centers (). Using these data, normalize the frequency counts using the overall area under the histogram. Plot this normalized histogram and overlay the theoretical PDF for the chosen parameters.

%For those who don't have access to 'histogram' function
%get un-normalized values from hist function with same number of bins as histogram function
numBins=50; %choose appropriately
[f,x]=hist(R,numBins); %use hist function and get unnormalized values
figure; plot(x,f/trapz(x,f),'b-*');%plot normalized histogram from the generated data

X = -4:0.1:4; %range of x to compute the theoretical pdf
fx_theory =   pdf('Normal',X,mu,sigma); %theoretical normal probability density
hold on; plot(X,fx_theory,'r'); %plot computed theoretical PDF
title('Probability Density Function'); xlabel('values - x'); ylabel('pdf - f(x)'); axis tight;
legend('simulated','theory');

Step 3: Theoretical PDF:

The given code snippets above,  already include the command to plot the theoretical PDF by using the ‘pdf’ function in Matlab. It you do not have access to this function, you could use the following equation for computing the theoretical PDF

The code snippet for that purpose is given next.

X = -4:0.1:4; %range of x to compute the theoretical pdf
fx_theory = 1/sqrt(2*pi*sigma^2)*exp(-0.5*(X-mu).^2./sigma^2);
plot(X,fx_theory,'k'); %plot computed theoretical PDF

Note:  The functions – ‘random’ and ‘pdf’ , requires statistics toolbox.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] John Mount, ‘Six Fundamental Methods to Generate a Random Variable’, January 20, 2012.↗
[2] Thomas, D. B., Luk. W., Leong, P. H. W., and Villasenor, J. D. 2007. Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages DOI = 10.1145/1287620.1287622 http://doi.acm.org/10.1145/1287620.1287622.↗

Topics in this chapter

Random Variables - Simulating Probabilistic Systems
● Introduction
Plotting the estimated PDF
● Univariate random variables
 □ Uniform random variable
 □ Bernoulli random variable
 □ Binomial random variable
 □ Exponential random variable
 □ Poisson process
 □ Gaussian random variable
 □ Chi-squared random variable
 □ Non-central Chi-Squared random variable
 □ Chi distributed random variable
 □ Rayleigh random variable
 □ Ricean random variable
 □ Nakagami-m distributed random variable
Central limit theorem - a demonstration
● Generating correlated random variables
 □ Generating two sequences of correlated random variables
 □ Generating multiple sequences of correlated random variables using Cholesky decomposition
Generating correlated Gaussian sequences
 □ Spectral factorization method
 □ Auto-Regressive (AR) model

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Linear Models – Least Squares Estimator (LSE)

Key focus: Understand step by step, the least squares estimator for parameter estimation. Hands-on example to fit a curve using least squares estimation

Background:

The various estimation concepts/techniques like Maximum Likelihood Estimation (MLE), Minimum Variance Unbiased Estimation (MVUE), Best Linear Unbiased Estimator (BLUE) – all falling under the umbrella of classical estimation – require assumptions/knowledge on second order statistics (covariance) before the estimation technique can be applied. Linear estimators, discussed here, do not require any statistical model to begin with. It only requires a signal model in linear form.

Linear models are ubiquitously used in various fields for studying the relationship between two or more variables. Linear models include regression analysis models, ANalysis Of VAriance (ANOVA) models, variance component models etc. Here, one variable is considered as a dependent (response) variable which can be expressed as a linear combination of one or more independent (explanatory) variables.

Studying the dependence between variables is fundamental to linear models. For applying the concepts to real application, following procedure is required

  1. Problem identification
  2. Model selection
  3. Statistical performance analysis
  4. Criticism of the model based on statistical analysis
  5. Conclusions and recommendations

Following text seeks to elaborate on linear models when applied to parameter estimation using Ordinary Least Squares (OLS).

Linear Regression Model

A regression model relates a dependent (response) variable y to a set of k independent explanatory variables {x1, x2 ,…, xk} using a function. When the relationship is not exact, an error term e is introduced.

If the function f is not a linear function, the above model is referred as Non-Linear Regression Model. If f is linear, equation (1) is expressed as linear combination of independent variables xk weighted by unknown vector parameters θ = {θ1, θ2,…, θk } that we wish to estimate.

Equation (2) is referred as Linear Regression model. When N such observations are made

where,
yi – response variable
xi – independent variables – known expressed as observed matrix X with rank k
θi – set of parameters to be estimated
e – disturbances/measurement errors – modeled as noise vector with PDF N(0, σ2 I)

It is convenient to express all the variables in matrix form when N observations are made.

Denoting equation (3) using (4),

Except for X which is a matrix, all other variables are column/row vectors.

Ordinary Least Squares Estimation (OLS)

In OLS – all errors are considered equal as opposed to Weighted Least Squares where some errors are considered significant than others.

If is a k ⨉ 1 vector of estimates of θ, then the estimated model can be written as

Thus the error vector e can be computed from the observed data matrix y and the estimated as

Here, the errors are assumed to be following multivariate normal distribution with zero mean and standard deviation σ2.

To determine the least squares estimator, we write the sum of squares of the residuals (as a function of ) as

The least squares estimator is obtained by minimizing . In order to get the estimate that gives the least square error, differentiate with respect to and equate to zero.

Thus, the least squared estimate of θ is given by

where the operator T denotes Hermitian Transpose (conjugate transpose).

Summary of computations

  1. Step 1: Choice of variables. Choose the variable to be explained (y) and the explanatory variables { x1, x2 ,…, xk } where x1 is often considered a constant (optional) that always takes the value 1 – this is to incorporate a DC component in the model.
  2. Step 2: Collect data. Collect n observations of y and for a set of known values of { x1, x2 ,…, xk }. Example: { x1, x2 ,…, xk } is the pilot data in OFDM using which we would like to estimate the channel impulse response θ and y is the received vector of samples. Store the observed data y in an – n⨉1 vector and the data on the explanatory variables in the n⨉k matrix X.
  3. Step 3: Compute the estimates. Compute the least squares estimates by the formula

The superscript T indicates Hermitian Transpose (conjugate transpose) operation.

Key Points

  • We do not need a probabilistic assumption but only a deterministic signal model.
  • It has a broader range of applications.
  • Least squares is unbiased.
  • Estimating the disturbance variance (k variables to estimate and n observations are available).
  • To keep the variance low, the number of observations must be greater than the number of variables to estimate.
  • The observation matrix X should have maximum rank – this leads to independent rows and columns which always happens with real data. This will make sure (XTX) is invertible.
  • Least Squares Estimator can be used in block processing mode with overlapping segments – similar to Welch’s method of PSD estimation.
  • Useful in time-frequency analysis.
  • Adaptive filters are utilized for non-stationary applications.

LSE applied to curve fitting

Matlab snippet for implementing Least Estimate to fit a curve is given below.

x = -5:.1:5; % set of x- values - known explanatory variables
y = 5.3 + 1.2* x; % Straight line without noise
e=randn(size(y));
y = y + e; % adding random noise to get observed variable - 
%Linear model - Y=Xa+e where a - parameters to be estimated

X = [ ones(length(x),1) x']; %first column treated aas all ones since x_1=1
y = y'; %column vector for proper dimension during multiplication
a = inv(X'*X)*X'*y  % Least Squares Estimator - equivalent code X\y
h=plot ( x , y , 'o'); %original data
hold on;
plot( x , a(1)+ a(2)*x , 'r-' ); %Fitted line
legend('observed samples',['y=' num2str(a(1)) '+' num2str(a(2)) 'x']) 
title('Least Squares Estimate for Curve Fitting');
xlabel('X values');
ylabel('Y values');

Simulation Results

Figure 1: Least Squares Estimate for Curve Fitting

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Related topics:

[1]An Introduction to Estimation Theory
[2]Bias of an Estimator
[3]Minimum Variance Unbiased Estimators (MVUE)
[4]Maximum Likelihood Estimation
[5]Maximum Likelihood Decoding
[6]Probability and Random Process
[7]Likelihood Function and Maximum Likelihood Estimation (MLE)
[8]Score, Fisher Information and Estimator Sensitivity
[9]Introduction to Cramer Rao Lower Bound (CRLB)
[10]Cramer Rao Lower Bound for Scalar Parameter Estimation
[11]Applying Cramer Rao Lower Bound (CRLB) to find a Minimum Variance Unbiased Estimator (MVUE)
[12]Efficient Estimators and CRLB
[13]Cramer Rao Lower Bound for Phase Estimation
[14]Normalized CRLB - an alternate form of CRLB and its relation to estimator sensitivity
[15]Cramer Rao Lower Bound (CRLB) for Vector Parameter Estimation
[16]The Mean Square Error – Why do we use it for estimation problems
[17]How to estimate unknown parameters using Ordinary Least Squares (OLS)
[18]Essential Preliminary Matrix Algebra for Signal Processing
[19]Why Cholesky Decomposition ? A sample case:
[20]Tests for Positive Definiteness of a Matrix
[21]Solving a Triangular Matrix using Forward & Backward Substitution
[22]Cholesky Factorization - Matlab and Python
[23]LTI system models for random signals – AR, MA and ARMA models
[24]Comparing AR and ARMA model - minimization of squared error
[25]Yule Walker Estimation
[26]AutoCorrelation (Correlogram) and persistence – Time series analysis
[27]Linear Models - Least Squares Estimator (LSE)
[28]Best Linear Unbiased Estimator (BLUE)

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Block Interleaver Design for RS codes

Introduction:

A Reed Solomon (RS) encoder, takes user data symbols and converts it into a n symbol wide codeword, by adding parity symbols. The error correcting capability of the RS code is computed as . That is, a RS code with parity symbols can correct a burst error of upto symbol errors.

Block Interleavers:

Suppose, assume that the dominant error mechanism in a channel is of burst type. A burst of length b is defined as a string of b unreliable consecutive symbols. If the expected burst length, b is less than or equal to t (the number of correctable symbol errors by RS coding), the code can be used as it is. However, if bursts length , the error correcting code will fail. This is where interleaving comes to our rescue.

Let us assume that and , where d (the interleaving depth) is an integer. The Reed-Solomon code can be used if we can spread the burst error sequence over several code blocks so that each block has no more than t errors (which can then be corrected). This can be accomplished using block interleaving as follows. Instead of encoding blocks of k symbols and then sending the encoded symbols consecutively, we can interleave the encoded blocks and transmit the interleaved data. In the case where ,the data bytes output from the Reed-Solomon encoder would appear as shown below , where bytes numbered 0 to 234 are the data bytes and bytes 235 to 254 are the parity check bytes.

Code Structure for Reed Solomon (RS) Codes

Here, the data is written row by row and read back column by column.Consider now the effect of a burst error of length , (where is the number of correctable errors per block) and for some , on the received symbols in the table. Because of the order in which the symbols are sent, a burst length less than or equal to will effect at most consecutive columns of the table, depending on where the burst starts. Notice that any single row (which corresponds to a codeword) has no more than v errors in it. If , these errors are within the error correction capability of the code and can be corrected. In this case, becomes the interleaving depth. The trade-off is that extra buffer space is required to store the interleaver table and additional delay is introduced. The worst case burst length determines the size of the table (and the interleaving depth) and the table size determines the amount of buffer space required and the delay.

Design Example:

Consider a (255,235) Reed Solomon coding system. This code can correct upto symbols errors. Lets assume that the channel that we are going to use, is expected to cause symbols. Then the interleaver depth is calculated as

In this case , an interleaver depth of 26 is enough to combat the burst errors introduced by the channel. The block interleaver dimensions would be (26 rows by 255 columns).

Matlab Code:

A sample matlab code that simulates the above mentioned block interleaver design is given below. The input data is a repeatitive stream of following symbols – “THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_“. A (255,235) Reed Solomon decoder (with correction capability of 10 symbols) is used. We assume that the channel is expected to produce a maximum of consecutive 20 symbols of burst error. The burst errors are denoted by ‘*‘.

%Demonstration of design of Block Interleaver for Reed Solomon Code
%Author : Mathuranathan for https://www.gaussianwaves.com
%License - Creative Commons - cc-by-nc-sa 3.0
clc;
clear;
%____________________________________________
%Input Parameters
%____________________________________________
%Interleaver Design for Reed-Solomon Codes
%RS code parameters
n=255;  %RS codeword length
k=235;  %Number of data symbols
b=20;  %Number of symbols that is expected to be corrupted by the channel
%____________________________________________

p=n-k;  %Number of parity symbols
t=p/2; %Error correction capability of RS code

fprintf('Given (%d,%d) Reed Solomon code can correct : %d symbols \n',n,k,fix(t));
fprintf('Given - expected burst error length from the channel : %d symbols \n',b);
disp('____________________________________________________________________________');
if(b>t)
    fprintf('Interleaving  MAY help in this scenario\n');
else
    fprintf('Interleaving will NOT help in this scenario\n');
end
disp('____________________________________________________________________________');
D=ceil(b/t)+1; %Intelever Depth

memory = zeros(D,n); %constructing block interleaver memory

data='THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_'; % A constant pattern used as a data

%If n>length(data) then repeat the pattern to construct a data of length 'n'
data=char([repmat(data,[1,fix(n/length(data))]),data(1:mod(n,length(data)))]);

%We sending D blocks of similar data
intlvrInput=repmat(data(1:n),[1 D]);

fprintf('Input Data to the Interleaver -> \n');
disp(char(intlvrInput));
disp('____________________________________________________________________________');

%INTERLEAVER
%Writing into the interleaver row-by-row
for index=1:D
    memory(index,1:end)=intlvrInput((index-1)*n+1:index*n);
end
intlvrOutput=zeros(1,D*n);
%Reading from the interleaver column-by-column
for index=1:n
    intlvrOutput((index-1)*D+1:index*D)=memory(:,index);
end

%Create b symbols error at 25th Symbol location for test in the interleaved output
%'*' means error in this case
intlvrOutput(1,25:24+b)=zeros(1,b)+42;
fprintf('\nInterleaver Output after being corrupted by %d symbol burst error - marked by ''*''->\n',b);
disp(char(intlvrOutput));
disp('____________________________________________________________________________');

%Deinteleaver
deintlvrOutput=zeros(1,D*n);
%Writing into the deinterleaver column-by-column
for index=1:n
    memory(:,index)=intlvrOutput((index-1)*D+1:index*D)';
end

%Reading from the deinterleaver row-by-row
for index=1:D
    deintlvrOnput((index-1)*n+1:index*n)=memory(index,1:end);
end
fprintf('Deinterleaver Output->\n');
disp(char(deintlvrOnput));
disp('____________________________________________________________________________');

Simulation Result:

Given : (255,235) Reed Solomon code can correct : 10 symbols
Given : expected burst error length from the channel : 20 symbols


Interleaving MAY help in this scenario


Input Data to the Interleaver ->
THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_T
HE_LAZY_DOG_
THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_
JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUIC
K_BROWN_FOX_JUMPS_OVER_THE_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_
QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_
LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JU
MPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUICK_BROWN_FOX
JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUIC
K_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY

DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_O
VER_THE_
____________________________________________________________________________

Interleaver Output after being corrupted by 20 symbols burst error – marked by ‘*‘->
TTTHHHEEE___QQQUUUIIICCC********************N___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVV
VEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK
___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHH
EEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWN
NN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___
DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJ
UUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEE
E___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOV
VVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___
BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE__
_
____________________________________________________________________________

Deinterleaver Output->
THE_QUIC*******_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_
LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUM
PS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BR
OWN_FOX_JUMPS_OVER_THE_THE_QUIC*******_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BRO
WN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_TH
E_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_
LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUIC******N_FOX_JUMPS_OVER_THE_L
AZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS
_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN
_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_
_______________
_____________________________________________________________

As we can see from the above simulation that, eventhough the channel introduces 20 symbols of consecutive burst error (which is beyond the correction capability of the RS decoder), the interleaver/deinterleaver operation has effectively distributed the errors and reduced the maximum burst length to 7 symbols (which is easier to correct by (255,235) Reed Solomon code.

See also:

[1] Introduction to Interleavers and deinterleavers
[2] Random Interleavers

Additional Resources:

[1] Notes on theory and construction of Reed Solomon Codes – Bernard Sklar
[2] Concatenation and Advanced Codes – Applications of interleavers- Stanford University