Modelling a Communication System and Bit Error Rate

In this example, we will see how we can model a simple communication channel. We will consider two cases : A simple Additive Gaussian noise channel and Rayleigh fading channel. Our signals are complex variables here which signifies both I and Q components.

Fig1: Communication System

Channel
A channel is the medium between the transmitter and receiver. The signals (for example electromagnetic waves) transmitted from the source face different materials in the medium (air, water, buildings, dust or static charges around) it may increase the signal strength or decrease it by some random factor. This is called fading and this factor is called fading variable.
It can also add some random component to the transmitted signal. We can see this as additive noise. When this component follows a Gaussian distribution we call it additive Gaussian Noise.

Signal Modelling
Let us generate a random bit stream which signifies the date need to be sent. We choose the probability of 1 occurring is 0.6. First we generate a random uniform random vector.

u=rand(1000,1)           

Now let us turn  all the elements bigger than 0.4 to 1 and rest to 0. Note that since the original vector is uniform between 0 and 1, the output vector will have 1 with 0.6 probability.

d=zeros(1000,1);
u(d>0.4)=1;

Modulation:

We can further modulate this signal using any modulation scheme. Let us choose few of modulation schemes and compare their outputs. 

a) BPSK

This means -1 denotes 0 bit and +1 denotes 1 bit. let us change all the 0s to -1 and 1s to +1. The transmitted signal can we written as

y=zeros(1000,1);
y(d==0)=-1;
y(d==1)=+1;

Note that each -1 and 1 denotes a sin pulse with phase -180 and +180 with time length equal to symbol duration . We will not consider the anolog signal which is in form concatenated sin signals. The receiver considers whole duration pulse as one and decide the phase. This requires synchronization which is not covered in this tutorial.

B) QPSK or QAM

Let us first divide the signal into I and Q data streams components. Each even numbered data bit is for Q.

dI=d(1:2:999);
dQ=d(2:2:1000);

The transmitted signal can we written as

x=zeros(500,1);x(dI==0 && dQ==0)=sqrt(1/2)*(1+1j);
x(dI==0 && dQ==1)=sqrt(1/2)*(-1+1j);
x(dI==1 && dQ==1)=sqrt(1/2)*(-1-1j);
 x(dI==1 && dQ==0)=sqrt(1/2)*(1-1j);

Channel Modeling

1. Additive gaussian Noise

Let us first consider the simple channel which just adds the Gaussian noise. noise vector is of the same dimension , one noise value for each symbol duration and each noise value  is independent of  previous
symbol. We choose the noise variance as 0.01 and mean as 0.

m=0;
s=sqrt(.01);
N=size(x);
n=m+s*rand(N);

So the output signal which is received at receiver is 

y=x+n;

2. Fading

Let us consider the fading which multiplies to the transmitted signal. As previous, here also each value is independent. We consider Rayleigh fading which means each symbol is Rayleigh distributed. .

h=sqrt(1/2)(randn(N,1)+1j*randn(N,1));
y=h*x+n; 

Receiver Model

We implement simple decoder which computes the phase of signal (as both of them are PSK) and decide upon this. For BPSK, it is equivalent to the following rule

dhat= 1 if y >0
dhat =0 if y<0

for QPSK, the rule is 

dIhat = 1 if imaginary part of y<0 otherwise dIhat=1;

dQhat = 1 if real part of y<0 otherwise dQhat=1;  

Let us implement both of these

BPSK:

dhat=zeros(1000,1);
dhat(y>0)=1;

QPSK:

dIhat=zeros(500,1);
dIhat(imag(y)<0)=1;
dQhat=zeros(500,1);
dQhat(real(y)<0)=1;
dhat=zeros(1000,1);
dhat(1:2:999)=dIhatl
dhat(2:2:1000)=dQhat;

BER performance

Now let us compare the bit error rate which is defined as the ratio of error bit to the total bits transmitted.

e=(d!=dhat);
ber=sum(e)/1000; 

Impulse Response: Is it just response to impulse

Hey All

After a long delay finally a post!! Today we will see why Impulse response is so important while rare cases are those where we will actually feed the impulse as input. We will take discreet LTI systems here for the understanding. While most of the things are true for continuous LTI systems too. There are some prerequisite for better understanding which I hope you all already know about. Here are they:
1. Discreet Time Fourier Transform \(X(e^{j\omega})\)
2. LTI - Linear Time Invariant Systems

Impulse Inputs and System's Response to it
Impulse are something which are just for a moment. Like You get hit by a ball while watching cricket in the stadium. Or the sudden breaks on the car. Fig 1 describes an impulse.  So Impulse is a signal which is one at a single time and zero everywhere else. We denote an Impulse by \(\delta[n]\) . It means that it is one only when \(n=0\). Similarily \(\delta[n-2]\) represent an impulse which is one only at \(n=2\).

Fig 1 Impulse Input

Although its response may not be impulse. You can feel the pain for a long time. or the car slows down, wiggles for a few seconds and then maintains its speed again. This Fig 2 shows a typical response. Let us call it  \(h[n]\).

Fig 2 Impulse Response: system's response to impulse

 Response (Output) of a system when an Impulse inputs is applied. Similarly a response to a step input is known as step response. Since system is LTI, it does not matter when or how you apply the impulse input. If it is given an impulse, it will always respond the same. and that is called impulse response which is a property of a system. But So is the step response. So why impulse response is treated in such an important fashion.

Response to a two rank input
Let us take an example with a response which is not impulse but looks like something shown in figure 3.

Fig 3 rank 2 input x[n]

$$x[0]=1$$ $$x[1]=0.4$$
What is the response to this input? Well we know the system is linear , so why not break the input to linear combinations of some input (basis inputs if you call) for which we know the input. If you think, you will find out impulse are such a choice which are simpler to find. ( Think in terms of finding the coefficients of a point in n dimensional in standard basis).
Fig 4 shows such a break up. Easy.

Fig 4 Breaking up the input into impulses

 So we can say $$x[n]=1\delta[n]+0.4\delta[n-1]$$
We can also write this as $$x[n]=x[0]\delta[n]+x[1]\delta[n-1]$$
Now use the linearity property.

\(\delta[n]\) results in \(h[n]\).
\(\delta[n-1]\) results in \(h[n-1]\).
So \(1\delta[n]+0.4\delta[n-1]\) results in \(1.h[n]+0.4h[n-1]\).
Fig 5 describes it.

Fig 5 Response to Each Impulse Add to Final Response

Woah!! So we can use impulse response even if the input is not impulse? Yes, the output is always an linear combination of impulse response.
Let us write the result in more generalized way
if \(x[n]\) is \(x[0]\delta[n]+x[1]\delta[n-1]\)
     Output (Response) is \(y[n]=x[0]h[n]+x[1]h[n-1]\)

Generalized Case:
What if the input have k+1 nonzero values at n=0,1,...,k. You can again write $$x[n]= x[0]\delta[n]+x[1]\delta[n-1]+...+x[k]\delta[n-k]$$
So the response is
$$y[n]=x[0]h[n]+x[1]h[n-1]+...+x[k]h[n-k]=\sum_{i=0}^{k}x[i]h[n-i]$$

So any response is a linear combination of shifted impulse responses with coefficients equal to signal values. That is why impulse response are so important and treated as basic property of a LTI system. If we know impulse response, we know  the response to any input.

Step Response
Let us try our findings on step input. For a step signal x[n] is always one for \(x\ge 0\). So response (or step response because it is response to step input) is $$y[n]=\sum_{i=0}^{\infty}h[n-i]$$. Sum (or integration) of h[n].

How to get back the impulse response if we know the response for an arbitrary input
Now this is little tricky. If you observe the equation $$y[n]=\sum_{i=0}^{k}x[i]h[n-i]$$ closely, you will know it is the convolution of x[n] and h[n].$$y[n]=x[n]*h[n]$$.

 So What do you do when you see convolution? Yes go to the fourier Domain. Let us denote DTFT of x[n],h[n] and y[n] as \(X(e^{j\omega}),H(e^{j\omega}) \) and \(Y(e^{j\omega})\). Here again \(X(e^{j\omega})\) is a system property because it is fixed and is known as Frequency response. So the earlier expression becomes $$Y(e^{j\omega})=X(e^{j\omega})H(e^{j\omega})$$ $$H(e^{j\omega})=\frac{Y(e^{j\omega})}{X(e^{j\omega})}$$. Remember that this function is constant regardless of what x[n] or y[n] you take if the system is LTI.

Book Listed on Mathworks

The book got reviewed by MathWorks Book Team and got listed on its website.

Do Checkout the online listing by yourself
http://www.mathworks.com/support/books/book62199.html

Other Online Listing Includes:
Authors Webpage:
https://sites.google.com/site/gkrabhishek/books

FaceBook Page:
http://www.facebook.com/pages/Matlab-By-Examples-Book/158144800875598

GooglePlus Page:
https://plus.google.com/109631828208021263344/posts

Median Filtering using MATLAB

In this blog we will learn how we can do a time domain filtering for an image. We will consider median filtering for this example. We will see two methods - first one is the iterative methhod which is time consumingg. Second methods uses matlab's vectorization and performs well.

Example Image
For this blog, we will take a very short image to understand how filtering is being done. Our sample image is shown in Fig1 (7x7).

Fig1 : Image A (magnified)

A =
8 5 8 0 6 8 7
9 9 1 8 1 6 7
1 9 4 9 7 3 1
9 1 9 6 0 9 4
6 9 7 7 2 0 4
0 9 9 7 0 4 6
2 4 6 3 0 3 7

You can create this by the following command

A=[8 5 8 0 6 8 7;9 9 1 8 1 6 7;1 9 4 9 7 3 1;9 1 9 6 0... 9 4; 6 9 7 7 2 0 4;0 9 9 7 0 4 6;2 4 6 3 0 3 7];

In actual cases you will read a image by using imread

A=imread('C:\users\Myname\codes\image1.jpg');
A=rgb2gray(A);

[l m]=size(A);

Kernel
We would be filtering with 3x3 kernel. It means that for each pixel of the output image, we would take a 3x3 neighborhood around the respective pixel in the original image, and assign the median color value in this neighborhood to the pixel of output image.

Output
Output would be a image of 7x7 too. Let us preallocate it with all zeros.

B= zeros(l,m);

Iterative method:
For each i and j, we first extract a 3x3 matrix around this point as shown in figure.

Fig2: Image Kernel Elaboration

K=A(i-1:i+1,j-1:j+1);

This will select a 3x3 matrix containing elements from i-1,i,i+1 row and j-1,j,j+1 column and assign it to K.

Now calculate median of this K and assign it to B(i,j);

B(i,j)=median(K);

We run a loop on i,j = 2 to 6. We cannot start with 1 because 3x3 kernel should fit inside the image to avoid any matrix index out of bound error. In the same way we cannot have i or j = 7. Note that we can take only those i where i-1, i,i+1 all three should come in the range 1 to 7. So with the j.

for i=2:6

  for j=2:6

    K= A(i-1:i+1,j-1:j+1);

    B(i,j)=median(K);

  end

end

 So for example for i=3,j=4, we take the kernel consisting of 2-3-4 row and 4-5-6 column as shown in figure 2. Elements are 1,4,6,7,7,9,9,9,9 and the median is 7 , so the output image's B(3,4)=7.

 This is done for all the indexes.

Padding:
We will first pad the original image to take care of the corner cases. To pad we will create a new matrix with l+2, m+2 dimension and assign the middle region from i=2 to l+1, j=2 to m+1 from the original image.

Y=zeros(l+2,m+2);
Y(2:l+1,2:m+1)=A;

Now remember rth,sth element of Y is equal to ith,jth element of A where i=r-1,j=s-1 for r,s>2.
Now you can go again i=2:l+1, 2:m+1, but remember to assign it as following

B(i-1,j-1)=median(K);

Vectorized Method:
The iterative method is little slow. MATLAB is optimized for vectorized operation. Let us try the same with vectorized way but it is little tricky. As you know, to calculate median, you need to pass MATLAB a matrix.

First we will create 9 shifted version of A (let us say A1, A2,... A9) so that for any i,j,A1(i,j) to A9(i,j) represent the 9 values in 3x3 neighborhood of A(i,j).
since we are considering i=2 to 6, these 9 matrixes will be 5x5 size;

A1=A(1:5,1:5);
A2=A(1:5,2:6);
A3=A(1:5,3:7);


A4=A(2:6,1:5);
A5=A(2:6,2:6);
A6=A(2:6,3:7);


A7=A(3:7,1:5);
A8=A(3:7,2:6);
A9=A(3:7,3:7);

Now we need to find median of 9 values A1(i,j) to A9(i,j) for all i,j. We first combine them into a 3D matrix and calculate the median in its first dimension.

C(1,:,:)=A1;
C(2,:,:)=A2;
C(3,:,:)=A3;
C(4,:,:)=A4;
C(5,:,:)=A5;
C(6,:,:)=A6;
C(7,:,:)=A7;
C(8,:,:)=A8;
C(9,:,:)=A9;
D=median(C,1);
B=squeeze(D);

And Done! Here is you answer B. This is pretty fast specially for large matrixes. There are only 9 operations (equal to kernel size) are happening, while in iterative one, m*n operations are happening, so if m=n=1000, there are 1000000 operations !!

Multi-variate Model Fitting using Taylor Series Method

In the previous example, we have learnt how to calculate a singlevariate model using a taylor series method when no knowledge about form of relation is known between input and output. This example will extend the same using iterative application of singlevariate model fitting.

We will consider a bivariate model and same can be easily generalized for higher dimensions.

Problem Statement:
We consider a system where there are two inputs/parameters (x and y) which we can change and observe the output (z). Our aim is to find a relation between this output and inputs but again we don't know what relation form(whether it is exponential, logarithmic, linear or trigonometric) exists.

Taylor Series Method
see this article for Taylor Series Method http://matlabbyexamples.blogspot.com/2011/12/singlevariate-model-fitting-for.html.

Iterative Taylor Series Method
We will first make one parameter (y) constant (let us say at y0) and change the other one (x) only. Then we try to fit a taylor series model in to it and we get coefficient b_n|at y=y0 such that  \[ f(x) = \sum_n{b_n x^n} \].
Now we will change the parameter y to y1 and do the same fit again to get a different set of coefficients  b_n|at y=y1. Do it for all value of y and each time you get a 1xN vector b containing b_n's. Now stack them such that  we get a matrix B each row represent b for one value of y. Now observe that ith column (let us call it B_i) represents b_i's for all values of y.

Now it is clear that these columns are only dependent on value of y. So again we can fit a taylor series to each of these columns such that  \[ b_n = \sum_m{c_m y^m} \].  This can be done by solving \[B_n = Y*C\]

So we have the final model as
\[ f(x) = \sum_n{\sum_m{c_{m,n} y^m} x^n} \].

MATLAB Implementation:

function C=taylorseriesfit2D(X,Y,Z);
%assumes X and Y are meshgrid format matrix
for i=1:length(y)
z=Z(i,:)';
B(i,:)=taylorseriesfit1D(X(1,:),z);
end

for i=1:size(B,2)
C(:,i)=taylorseriesfit1D(Y(:,1),B(:,i));
end

where taylorseriesfit1D is given as

function B= taylorseriesfit1D(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7];
B=Xmat\z;

Simulation Data generation:

x=[0:0.1:1]';
y=[0:.1:pi/2]';
[X Y ]=meshgrid(x,y);
Z=exp(X).*exp(Y);
C=taylorseriesfit2D(X,Y,Z)

C will come like for N=7.

1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
0.5000 0.5000 0.2500 0.0833 0.0208 0.0043 0.0006 0.0002
0.1669 0.1669 0.0835 0.0278 0.0069 0.0014 0.0002 0.0001
0.0410 0.0410 0.0205 0.0068 0.0017 0.0004 0.0000 0.0000
0.0092 0.0092 0.0046 0.0015 0.0004 0.0001 0.0000 0.0000
0.0008 0.0008 0.0004 0.0001 0.0000 0.0000 0.0000 0.0000
0.0004 0.0004 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000

Final Guesses 
Again we can guess the model by observing these coefficients and plotting them with respect to m,n (Fig1). First thing we can observer that it is symmetric, therefor exihibit the same relations for both the inputs to ouput.

Fig1: Plot C vrs m,n

Again, we cannot guess, even then these matrixes are sufficient for all numerical purposes.

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

In this example we will talk about how we can fit a model/function in a single variate simulation data. In many practical cases, we get some observation data after changing an input variable and our aim is to find the relation between input and output. If we know the relation form atleast (like y=mx+c or y=d*sin(x)^2), then we can find the values of these coeefficients y,c or d by linear/non linear regression. But what if we dont know the relation form, we dont have any information about whether they follow exponential form or logarthim relation. Here is a trick which we are going to discuss in this example.

Taylor Series fit
Let i/o relation be represented by a function z=f(x). We assume here that f(x) is continuous and differentiable, which is generally the case with practical/real time systems. Now each such function can be represented using a taylor's series
\[ f(x) = \sum_n{b_n x^n} \] we will limit the summation to n=N
Since we know the data x and z=f(x), we can calculate the coefficients bi .
Let \[B=\left[\begin{array}{ccccc} b_0 & b_1 & b_2  & ... & b_N \end{array} \right]^T \] and \[ Xmat =\left[ \begin{array}{ccccc} 1 & x^1 & x^2 & ... & x^N \end{array} \right] \]
so \[ z=f(x)= Xmat B \]

This is also true when x is a column vector mx1 and contains m data points and z contains values of f at those m point.
So B can be calculated as \[B=Xmat^{-1} Z\]

function: taylorseriesfit()
Writing the function in matlab is very easy. we just need to follow the steps as explained above.
Since our function need to matrixs x and z, these will be input to the function and B will be the output of the function. We will take N=4 here. Increasing the value may result in accuracy of your fit.

function B= taylorseriesfit(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x^4];
B=Xmat\z; 

Here x and z are column vectors representing data points.

Generation of Simulation Data 
We will generate a sample Simulation data for our example. In real example, you will get this data from some experiment after changing an input parameter(x) and observing the output(z).

x=[0:0.1:1]';

z=exp(x);

Calling the defined function

B=taylorseriesfit(x,z);
disp(B);

You will get B as

B =
1.0000
0.9988
0.5099
0.1400
0.0696

Now your model is
\[ z= x*B = B(1)+B(2)x+B(3)x^2+B(4)x^3+B(5)x^4 \] \[ z= 1+0.99x+0.50 x^2+0.14x^3+.07x^4 \]

Error
we can calculate the modelling error (MMSE) by comparing the actual output and the simulation output.

z_o = Xmat*B;
e=sqrt(sum((z-z_o).^2));
disp(e);

you will get something like

e= 6.9400e-005

which is pretty good.

Final Observation
Now if we observe the coefficients and plot them with respect with n, we may guess the actual function. Fig1 represents such a graph for N=7 case.

Fig1 : b_i plotted versus i

like for this example, it looks like 1/n!, so this is an exponential relations. Well, if we can guess, it is good. But if we cannot, still we have a relation which is sufficient for all numerical purposes.

Circuit Solver using MATLAB Programming

In this blog we will learn to simulate a circuit using MATLAB Programming which is the basic step to create your own circuit solver. The circuit we are going to solve is shown in Fig 1. Here V is 6V.

Fig 1: Circuit

Create NetList and SourceList
You need to first represent this circuit in terms of matrixes so that you can use them in MATLAB. These two matrixs will describe the network- NetList and SourceList
Let us first give a number to all the nodes in the circuit. Two nodes should be separated by atleast an element.
The considered circuit has 4 nodes and let us choose node 4 as the reference/ground node.

Now create a matrix named as NetList with three columns N1, N2 and R. Each row of this matrix will represent a R resistance element between Node N1 and N2.

Now create a matrix named SourceList with two columns N and V. Each row of this matrix will represent a N node with V voltage. Since only two node voltages are know, this will contain two rows.

Fig 2: NetList and Source List

Create A and B matrix
If Total Nodes are N
A is of N x N
Aij=
        Case 1: if i is not a Source Nodes
                       -1/Rij where Rij is the resistance between i and j (if i!=j)
                       Sum of inverse of all resistance connected to i (if i==j)

        Case 2: if i is a source Node
                        Aii=1 and
                        rest Aik=0

If Total Nodes are N
B is of N x 1
Bi=
       Case 1: if i is not a Source Node
                    0
       Case 2: if i is a source Node
                    Voltage of that source Node



Nodal Equations

Using the above nodal equation can be represented as
\[\mathbf{AV=B}\]\[\mathbf{V}=[V_1 V_2 ... V_N]^T\]\[V=A^{-1}B\]

This V represent the Voltage of all the nodes, thus circuit is solved.


MATLAB Code
Initialization of the matrixs

NetList=[1 2 2; 1 3 4; 2 3 5.2; 3 4 6; 2 4 3];
Vnod=[1 6; 4 2];
l=size(NetList,1);
N=max([NetList(:,1) ;
NetList(:,2)]);
A=zeros(N,N);
 B=zeros(N,1);

A Matrix Creation

for i=1:l
       n1=NetList(i,1);
       n2=NetList(i,2);
       if n1==n2
       else
          A(n1,n2)=A(n1,n2)-1/NetList(i,3);
          A(n2,n1)=A(n2,n1)-1/NetList(i,3);
          A(n1,n1)=A(n1,n1)+1/NetList(i,3);
          A(n2,n2)=A(n2,n2)+1/NetList(i,3);
      end
 end

B Matrix Creation

for i=1:size(Vnod,1)

        A(Vnod(i,1),:)=zeros(1,N);

        A(Vnod(i,1),Vnod(i,1))=1;

        B(Vnod(i,1),1)=Vnod(i,2);

end

Solution

Vo=A\B; 

disp(Vo);