# How to generate AWGN noise in Matlab/Octave (without using in-built awgn function)

(5 votes, average: 5.00 out of 5)

Looking for the proper way to generate AWGN noise in Matlab/Octave ? Here you go…

## AWGN – the in-built function

Matlab/Octave communication toolbox has an inbuilt function named – $$awgn()$$ with which one can add an Additive Gaussian White Noise to obtain the desired $$Signal-to-Noise Ratio$$ (SNR). The main usage of this function is to add AWGN to a clean signal (infinite SNR) in order to get a resultant signal with a given SNR (usually specified in dB). This usage is often found in signal processing/digital communication applications. For example, in Monte Carlo simulations involving modeling of modulation/demodulation systems, the modulated symbols at the transmitter are added with a random noise of specific strength, in order to simulate a specific $$E_b/N_0$$ or $$E_s/N_0$$ point.

The function $$y=awgn(x,SNR,’measured’)$$, first measures the power of the signal vector $$x$$ and then adds white Gaussian Noise to $$x$$ for the given $$SNR$$ level in dB. The resulting signal y is guaranteed to have the specified SNR.

## Custom function to add AWGN noise

If you do not have the communication toolbox, or if you would like to mimic the in-built AWGN function in any programming language, the following procedure can be used.
1) Assume, you have a vector $$x$$ to which an AWGN noise needs to be added for a given $$SNR$$ (specified in dB).
2) Measure the power in the vector $$x$$
$$E_s=\frac{1}{L} \sum_{i=0}^{L-1} \left|x[i] \right|^2 ; \;\;\; where \;\;\;\; L=length(x)$$
3) Convert given $$SNR$$ in dB to linear scale ($$SNR_{lin}$$) and find the noise vector (from Gaussian distribution of specific noise variance)  using the equations below

$$noise= \begin{cases} \sqrt{\frac{E_s}{SNR_{lin}}}*randn(1,L) & if \;\; x \;\; is \;\; real \\ \\ \sqrt{\frac{E_s}{2* SNR_{lin}}}*\left [randn(1,L)+j*randn(1,L) \right] & if \;\; x \;\; is \;\; complex \end{cases}$$

4) Finally add the generated noise vector to the signal $$x$$

$$y=x+noise$$

## The custom function – add_awgn_noise():

The custom function written in Matlab, that mimics the $$awgn$$ function is given below. It can be easily ported to Octave.

## Comparison & Testing:

Let’s cross check the results obtained from the above function with that of the standard in-built $$awgn$$ function in Matlab. Testing and comparison is done using two test waveforms – 1) sawtooth waveform (represented by a vector containing only real numbers) , 2) A complex sinusoidal waveform (vector having both real and imaginary part). The testing below is done for SNR=5 dB. Users are advised to test the function for various ranges of SNRs.

Test code and results for sawtooth waveform (the input signal x is a vector of real numbers) is given next. Test of linearity indicates that the results from the custom function matches with that of the in-built function.

Test code and results for complex sinusoidal waveform (the input signal x is a vector of complex numbers) is given next. Test of linearity indicates that the results from the custom function matches with that of the in-built function.

How to generate AWGN noise in Matlab/Octave by Mathuranathan Viswanathan is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.

# Power Delay Profile

(7 votes, average: 3.86 out of 5)

The Power Delay Profile $$p(\tau)$$ gives the distribution of signal power received over a multipath channel as a function of propagation delays.  It is obtained as the spatial average of the complex baseband channel impulse response as

$$p(\tau) = R_{hh}(0,\tau)= E\left [ \left | h(t,\tau) \right | ^2 \right ] \;\;\;\;\;\;\;\;\;\;\;\; (1)$$

As discussed in the previous post, it can also be derived from scattering function as given below

$$p(\tau) = \int_{-\infty }^{\infty } S(f,\tau) df \;\;\;\;\;\;\;\;\;\;\;\; (2)$$

where $$f$$ denotes Doppler Frequency, $$\tau$$ denotes multipath propagation delay, $$S(f,\tau)$$ denotes the scattering function.

In a Power Delay Profile plot, the signal power of each multipath is plotted against their respective propagation delays. A sample power delay profile plot, shown below, indicates how a transmitted pulse gets received at the receiver with different signal strength as it travels through a multipath channel with different propagation delays ($$\tau_0,\tau_1$$ and $$\tau_2)$$.

Power Delay Profile is usually supplied as a table of values obtained from empirical data and it serves as a guidance to system design. Nevertheless, it is not an accurate representation of the real environment in which the mobile is destined to operate at. For example, the 3GPP spec specifies the power delay profile for various environments as follows.

## Maximum Excess Delay Definition and applications

With Power Delay Profile, one can classify a multipath channel into frequency selective or frequency non-selective category. The derived parameter, namely, Maximum Excess Delay together with the symbol time of each transmitted symbol, can be used to classify the channel into frequency selective or non-selective channel.

Power Delay Profile can be used to estimated the average power of a multipath channel, measured from the first signal that strikes the receiver to the last signal whose power level is above certain threshold. This threshold is chosen based on receiver design specification and is dependent on receiver sensitivity and noise floor at the receiver.

Maximum Excess Delay, also called Maximum Delay Spread, denoted as $$T_m$$, is the relative time difference between the first signal component arriving at the receiver to the last component whose power level is above some threshold. Maximum Delay Spread $$T_m$$ and the symbol time peroid $$T_{sym}$$ can be used to classify a channel into frequency selective or non-selective category. This classification can also be done using Coherence Bandwidth (a derived parameter from Spaced Frequency Correlation Function which in turn is the frequency domain representation of power delay profile).

Maximum Excess Delay is also an important parameter in mobile positioning algorithm. The accuracy of such algorithm depends on how well the Maximum Excess Delay parameter conforms with measurement results from actual environment.

When a mobile channel is modeled as a FIR filter (tapped delay line implementation), as in CODIT channel model[1], the number of taps of the FIR filter is determined by the product of maximum excess delay and the system sampling rate.

The Cyclic Prefix in a OFDM system is typically determined by the maximum excess delay or by the RMS delay spread of that environment [2].

## Classification of Channel – signal spreading in Time domain:

A channel is classified as Frequency Selective, if the maximum excess delay is greater than the symbol time period, i.e, $$T_m > T_{sym}$$. This introduces Inter Symbol Interference (ISI) into the signal that is being transmitted, thereby distorting it.  This occurs since the signal components (whose power are above the threshold or the maximum excess delay)  due to multipath extend beyond the symbol time. ISI can be mitigated at the receiver by an equalizer.

In a frequency selective channel , the channel output can be expressed as the convolution of input signal and the channel impulse response plus some noise.

$$y(t) = h(t,\tau) \circledast x(t) + w(t) \;\;\;\;\;\;\;\;\;\;\;\; (3)$$

On the other hand, if the maximum excess delay is less than the symbol time period, i.e, $$T_m < T_{sym}$$, the channel is classified as frequency non-selective or “flat” channel. Here, all the scattered signal components (whose power are above the specified threshold or the maximum excess delay) due to the multipath, arrive at the receiver within the symbol time. This will not introduce any ISI, but the received signal is distorted due to inherent channel effects like SNR condition. Equalizers in the receiver are not needed.  A time varying non-frequency selective channel is obtained by assuming that the impulse response $$h(t,\tau) = h(t) \delta(\tau)$$. Thus the output of the channel can be expressed as

\begin{align*} y(t) &= h(t,\tau) \circledast x(t) + w(t) \\ &=\int h(t,\tau) x(t-\tau) d \tau + w(t) \\ &=\int h(t) \delta(\tau) x(t-\tau) d \tau + w(t) \\ &=h(t)x(t) + w(t) \end{align*} \;\;\;\;\;\;\;\;\;\;\;\; (4)

Note, that the output of the channel can be expressed simply as product of time varying channel response and the input signal. If the channel impulse response is a deterministic constant, i.e, time in-varying, then the non-frequency selective channel is expressed as follows by assuming $$h(t,\tau) = h \delta(\tau)$$

$$y(t) = h x(t) + w(t) \;\;\;\;\;\;\;\;\;\;\;\; (5)$$

This is the simplest situation that can occur. In addition to that, if the noise in the above equation is white Gaussian noise, the channel is called Additive White Gaussian Noise (AWGN) channel.

## Characterization of Frequency Selective Channels:

Average delay and the RMS delay spread are two most important parameters that characterize a frequency selective channel. They are derived from Power Delay Profile.

### Average delay:

Simply the statistical mean of the delay that a signal undergoes when transmitted over a multipath channel. For a frequency selective WSSUS channel, the average delay is equal to the first moment of the power delay profile $$p(\tau)$$. For a discrete channel it is calculated as

$$\bar{\tau}=\frac{\sum \tau p(\tau)}{\sum p(\tau)} \;\;\;\;\;\;\;\;\;\;\;\; (6)$$

If the given PDP values are continuous in terms of time delays (as given in the sample plot below), replace the summation with integral and integrate it with respect to $$d \tau$$.

### RMS delay spread:

RMS delay spread is equal to the second central moment of power delay profile $$p(\tau)$$. It is similar to the standard deviation of a statistical distribution. For a discrete channel it is given by

$$\sigma_{\tau}=\sqrt{\bar{\tau^2}-(\bar{\tau})^2} \;\;\;\;\;\;\;\;\;\;\;\; (7)$$

where,

$$\bar{\tau^2}=\frac{\sum \tau^2 p(\tau)}{\sum p(\tau)} \;\;\;\;\;\;\;\;\;\;\;\; (8)$$

If the given PDP values are continuous in terms of time delays (as given in the sample plot below), replace the summation with integral and integrate it with respect to $$d \tau$$.

The ratio of RMS delay spread and symbol time duration quantifies the strength of Inter Symbol Interference. This ratio determines the complexity of the equalizer required at the receiver. Typically, when the symbol time period is greater than 10 times the RMS delay spread, no ISI equalizer is needed in the receiver.

A sample power delay profile is plotted below. Average delay, RMS delay spread, and the maximum excess delay ( given a threshold of 100 dBm ) are all marked.

## References:

 [1] Andermo, P.G.; Larsson, G., “Code division testbed, CODIT,” Universal Personal Communications, 1993. Personal Communications: Gateway to the 21st Century. Conference Record., 2nd International Conference on , vol.1, no., pp.397,401 vol.1, 12-15 Oct 1993 [2] Huseyin Arslan, Cognitive Radio, Software Deﬁned Radio, and Adaptive Wireless Systems, pp. 238, 2007, Dordrecht, Netherlands, Springer.

# Simulation of M-PSK modulation techniques in AWGN channel

(2 votes, average: 3.50 out of 5)

A generic simulation technique to simulate all M-PSK modulation techniques (for upto $$M=32$$) is given here. The given simulation code is very generic, and it plots both simulated and theoretical symbol error rates for all M-PSK modulation techniques (upto $$M=32$$).

### M-PSK Modulation and simulation methodology:

The general expression for a M-PSK signal set is given by

$\dpi{120} {\color{DarkBlue}s_i \left( t \right)=V cos \left[ 2 \pi f_c t \; – \; \frac{ \left( i - 1 \right) 2 \pi }{M}\right] \;\;\; \text{ where i = 1,2,…,M} \;\; \rightarrow (1) }$

Here M defines the number of constellation points in the constellation diagram and essentially the M-PSK type. For example $$M=4$$ implies 4-PSK or QPSK, $$M=8$$ implies 8-PSK. The value of M depends on another parameter $$k$$ – the number of bits we wish to squeeze into a single M-PSK symbol. For example if we wish to squeeze in 3 bits ($$k=3$$) in one transmit symbol, $$M = 2^k = 2^3 = 8$$. This gives us 8-PSK configuration.

Equation (1) can be separated into cos and sin terms as follows

$\dpi{120} {\color{DarkBlue}s_i(t) = V cos \left[ \frac{(i-1)2\pi}{M} \right] cos(2 \pi f_c t) \; + \; V sin \left[ \frac{(i-1)2\pi}{M} \right] sin(2 \pi f_c t) \rightarrow \;\;\;(2) }$

This can be written as

$$\begin{matrix} s_i(t) = s_{i1}\phi_{1}(t) + s_{i2}\phi_{2}(t) \;\; \rightarrow (3)\\ where,\\ \\ s_{i1} = \sqrt{E_s}\; cos \left[ \frac{(i-1)2\pi}{M} \right] \;,\; s_{i2} = \sqrt{E_s}\; sin \left[ \frac{(i-1)2\pi}{M} \right]\\ \\ \phi_{1}(t) = \frac{V cos(2\pi f_c t)}{\sqrt{E_s}},\phi_{2}(t) = \frac{V sin(2\pi f_c t)}{\sqrt{E_s}}\end{matrix}$$

Here $$\phi_{1}(t)$$ and $$\phi_{2}(t)$$ are orthonormal basis functions that follows from Gram-Schmidt orthogonalization procedure and $$s_{i1}$$ and $$s_{i2}$$ are the coefficients of each signaling point in the M-PSK constellation. $$E_s$$ is the symbol energy usually normalized to $$1/\sqrt{2}$$ The constellation points on the M-PSK constellation lie $$2 \pi /M \;$$ radians apart and are placed on a circle of radius $$\sqrt{E_s}$$. The coefficients $$s_{i1}$$ and $$s_{i2}$$ are termed as inphase (I) and quadrature-phase (Q) components respectively.

The ideal constellation diagram for M-PSK contains M equally spaced signaling points that are located at the distance $$\sqrt{E_s}$$ from the origin. Figure 1 illustrates the ideal constellation diagram for 8-PSK constellation.

The generated I and Q components are then added with AWGN noise of required variance depending on the required Es/N0. The received signal’s constellation is corrupted with noise and the detection is based on comparing the received symbols with the ideal signaling points and making a decision based on the minimum distance.

Finally the simulated and theoretical symbol error rates are computed.

The theoretical symbol error rate [1] for M-PSK modulation is given by

\dpi{120} {\color{DarkBlue} \begin{align*} P_s &= 2Q\left[\sqrt{2E_s} \; sin\left(\frac{\pi}{M}\right)\right] \\ &= erfc\left [ \sqrt{E_s} \; sin\left(\frac{\pi}{M}\right)\right ] \end{align*} \;\;\;\;\; \rightarrow(4) }

### Matlab Simulation:

#### File 1: simulateMPSK.m

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

# Intuitive derivation of Performance of an optimum BPSK receiver in AWGN channel

(4 votes, average: 5.00 out of 5)

BPSK modulation is the simplest of all the M-PSK techniques. An insight into the derivation of error rate performance of an optimum BPSK receiver is essential as it serves as a stepping stone to understand the derivation for other comparatively complex techniques like QPSK,8-PSK etc..

Understanding the concept of Q function and error function is a pre-requisite for this section of article.

The ideal constellation diagram of a BPSK transmission (Figure 1) contains two constellation points located equidistant from the origin. Each constellation point is located at a distance $latex \sqrt{Es}$ from the origin, where Es is the BPSK symbol energy. Since the number of bits in a BPSK symbol is always one, the notations – symbol energy (Es) and bit energy (Eb) can be used interchangeably (Es=Eb).

Assume that the BPSK symbols are transmitted through an AWGN channel characterized by variance = N0/2 Watts. When 0 is transmitted, the received symbol is represented by a Gaussian random variable ‘r’ with mean=S0 = $latex \sqrt{Es}$ and variance =N0/2. When 1 is transmitted, the received symbol is represented by a Gaussian random variable – r with mean=S1= $latex \sqrt{Es}$ and variance =N0/2. Hence the conditional density function of the BPSK symbol (Figure 2) is given by,

 Figure 1: BPSK – ideal constellation Figure 2: PDF of BPSK Symbols

$latex \begin{matrix} f(r\mid 0_T)=\frac{1}{\sqrt{\pi N_0}}exp\left \{ -\frac{\left( \hat{r}-\hat{s}_0\right )^2}{N_0} \right \} \;\;\;\;\;\rightarrow (1A) \\ \\ f(r\mid 1_T)=\frac{1}{\sqrt{\pi N_0}}exp\left \{ -\frac{\left( \hat{r}-\hat{s}_1\right )^2}{N_0} \right \}\;\;\;\;\;\rightarrow (1B) \end{matrix} &s=1&fg=0000A0$

An optimum receiver for BPSK can be implemented using a correlation receiver or a matched filter receiver (Figure 3). Both these forms of implementations contain a decision making block that decides upon the bit/symbol that was transmitted based on the observed bits/symbols at its input.

When the BPSK symbols are transmitted over an AWGN channel, the symbols appears smeared/distorted in the constellation depending on the SNR condition of the channel. A matched filter or that was previously used to construct the BPSK symbols at the transmitter. This process of projection is illustrated in Figure 4. Since the assumed channel is of Gaussian nature, the continuous density function of the projected bits will follow a Gaussian distribution. This is illustrated in Figure 5.

After the signal points are projected on the basis function axis, a decision maker/comparator acts on those projected bits and decides on the fate of those bits based on the threshold set. For a BPSK receiver, if the a-prior probabilities of transmitted 0’s and 1’s are equal (P=0.5), then the decision boundary or threshold will pass through the origin. If the a-priori probabilities are not equal, then the optimum threshold boundary will shift away from the origin.

Considering a binary symmetric channel, where the a-priori probabilities of 0’s and 1’s are equal, the decision threshold can be conveniently set to T=0. The comparator, decides whether the projected symbols are falling in region A or region B (see Figure 4). If the symbols fall in region A, then it will decide that 1 was transmitted. It they fall in region B, the decision will be in favor of ‘0’.

For deriving the performance of the receiver, the decision process made by the comparator is applied to the underlying distribution model (Figure 5). The symbols projected on the $latex \phi$ axis will follow a Gaussian distribution. The threshold for decision is set to T=0. A received bit is in error, if the transmitted bit is ‘0’ & the decision output is ‘1’ and if the transmitted bit is ‘1’ & the decision output is ‘0’.

This is expressed in terms of probability of error as,

$\dpi{120} {\color{DarkBlue}P(error) = P( 1\; decided , 0\; transmitted )+P(0\; decided,1 \;transmitted) \rightarrow (2)}$
Or equivalently,

$latex P(error) = P(1_D,0_T) + P(0_D,1_T) \;\;\;\;\;\;\;\;\rightarrow (3) &s=1&fg=0000A0$

By applying Bayes Theorem, the above equation is expressed in terms of conditional probabilities as given below,

$latex P(error) = P(0_T) P( 1_D \mid 0_T ) + P(1_T) P(0_D \mid 1_T) \;\;\;\rightarrow (4) &s=1&fg=0000A0$

$\dpi{120} {\color{DarkBlue} P(error) = P(0_T)\int_{0}^{\infty}f(\hat{r}\mid0_T)d\hat{r} + P(1_T)\int_{-\infty}^{0}f(\hat{r}\mid1_T)d\hat{r} \rightarrow (5)}$
Since a-prior probabilities are equal P(0T)= P(1T) =0.5, the equation can be re-written as

$\dpi{120} {\color{DarkBlue}P(error) = 0.5\int_{0}^{\infty}f(\hat{r}\mid 0_T)d\hat{r} + 0.5\int_{-\infty}^{0}f(\hat{r}\mid 1_T)d\hat{r} \;\;\;\;\;\;\;\;\rightarrow (6)}$

Intuitively, the integrals represent the area of shaded curves as shown in Figure 6. From the previous article, we know that the area of the shaded region is given by Q function.

For the continuous density functions given in equation 1A, the mean=S0 = $latex \sqrt{Es}$  and variance =N0/2, the area of region of error in Figure 6a is given by

$latex P( 1_D \mid 0_T ) = Q\left (\sqrt{\frac{E_s}{N_0/2}}\right) \;\;\;\;\;\;\;\;\rightarrow (7) &s=1&fg=0000A0$

Similarly,

$latex P( 0_D \mid 1_T ) = Q\left(\sqrt{\frac{E_s}{N_0/2}}\right) \;\;\;\;\;\;\;\;\rightarrow (8) &s=1&fg=0000A0$

From (4), (6), (7) and (8),

\dpi{120} {\color{DarkBlue}\begin{align*} P(error) & = 0.5 Q\left(\sqrt{\frac{E_s}{N_0/2}}\right) + 0.5 Q\left(\sqrt{\frac{E_s}{N_0/2}}\right) \\ & = Q\left(\sqrt{\frac{2E_s}{N_0}}\right) = 0.5 \; erfc \left(\sqrt{\frac{E_s}{N_0}}\right) \end{align*} \;\;\rightarrow (9)}
For BPSK, since Es=Eb, the probability of symbol error (Ps) and the probability of bit error (Pb) are same. Therefore, expressing the Ps and Pb in terms of Q function and also in terms of complementary error function :

$latex P_s=P_b= Q\left(\sqrt{\frac{2E_s}{N_0}}\right) = Q\left(\sqrt{\frac{2E_b}{N_0}}\right) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\rightarrow (11) &s=1&fg=0000A0$
$latex P_s=P_b= 0.5 \; erfc\left(\sqrt{\frac{E_s}{N_0}}\right) = 0.5 \; erfc\left(\sqrt{\frac{E_b}{N_0}}\right) \;\;\rightarrow (12)&s=1&fg=0000A0$

# Eb/N0 Vs BER for BPSK over Rician Fading Channel

(4 votes, average: 3.00 out of 5)
This post is a part of the ebook : Simulation of digital communication systems using Matlab – available in both PDF and EPUB format.
In one of the previous article, simulation of BPSK over Rayleigh Fading channel was discussed. This article deals with simulation of another important fading type : Rician Fading.

Rayleigh Fading model is used to simulate environments that has multiple scatterers and no Line Of Sight (LOS) path. If there are sufficient multiple scatterers in the environment, all the reflected signals that appear at the receiver front end becomes uncorrelated in amplitude (this brings mean =0) and phase evenly distributed between 0 to 2π. Eventually the I (inphase) and Q (quadrature) phase components become Gaussian distributed (by central limit theorem) and their envelope $$Z = \sqrt{I^{2} + Q^{2}}$$ becomes Rayleigh distributed.
Note that distribution I2 is called Chi-Square distribution. If I2 has mean=0, then it is called central-Chi-square distribution and if mean ≠ 0 it is called non-central-Chi-square distribution

### Rayleigh Distribution:

Consider two Gaussian random variables X and Y with zero mean and equal variance σ2. Then the transformation $$Z = \sqrt{X^{2} + Y^{2}}$$ is Rayleigh Distributed and Z2 is exponentially distributed.

Rician Fading model is used to simulate environments that produces multiplath components and also a dominant Line Of Sight (LOS) component. The LOS component is called “specular” component and the multipath component is called “random or scatter” component. The amplitude distribution of the specular component will have non-zero mean, where as, the random component will have zero-mean.

### Rician Distribution:

Consider two Gaussian random variables X and Y. Here X models the specular component (LOS) and Y models the random/scatter component.By definition, X has non-zero mean (m), Y has zero mean and both have equal variance $$\sigma^2$$. Then the transformation, $$Z = \sqrt{X^{2} + Y^{2}}$$ is Rician Distributed.

The ratio of power of specular component to the power of random component is called Rician $$\kappa$$ factor and it is defined as

$$\kappa= \frac{m^2}{2 \sigma^{2}}$$

It can be immediately ascertained that Rayleigh Fading is related to central Chi Square distribution (due to zero mean) and the Rician Fading is related to non-central Chi Square distribution (due to non zero mean).

### Performance simulation of BPSK over Rician Fading Channel with AWGN noise:

The simulation strategy is very similar to the technique used for simulating Rayleigh Fading. The only difference here is that in Rician Fading, $$\kappa$$ factor is one of inputs that defines the relative strength of LOS component and the multiplath scattered component.

To simulate a Rician Fading channel, mean and sigma has to be calculated with the given Rician$$\kappa$$ factor.

Lets define mean ($$m$$ ) and sigma ($$\sigma$$) as given below, so that it satisfies the equation given above.

$$m = \sqrt{\frac{\kappa }{(\kappa +1)}}$$

$$\sigma = \sqrt{\frac{1}{2*(\kappa +1)}}$$

As discussed above, in Rician Fading, the specular component ($$X$$) has to be a Gaussian random variable with mean=$$m$$ and sigma=$$\sigma$$, but the scatter component ($$Y$$) has to be generated with mean=0 and sigma=$$\sigma$$.

In Matlab randn function generates Gaussian random numbers with mean=0 and sigma = 1. To generate $$X$$ component with mean=$$m$$ and sigma=$$\sigma$$, the output of randn has to be multiplied with $$\sigma$$ and added with m. To generate $$Y$$ component with sigma=σ, the output of the randn function has to be multiplied with σ. The rest of the simulation code is similar to that of BPSK over Rayleigh/AWGN channel simulation.

### Theoretical BER :

Theoretical BER for BPSK over Rician Fading Channel with AWGN noise is given by the following expression[1]

$$P_b = \frac{1}{2} erfc \left[ \sqrt{ \frac{ \kappa \left( E_b/N_0 \right) }{ \left( \kappa + E_b/N_0 \right) } } \right]$$

### Simulation Model:

The following model is used for the simulation of BPSK over Rician Fading channel (with AWGN noise).

### Matlab Code:

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

# Eb/N0 Vs BER for BPSK over Rayleigh Channel and AWGN Channel

(7 votes, average: 4.43 out of 5)

This post is a part of the ebook : Simulation of digital communication systems using Matlab – available in both PDF and EPUB format.

The phenomenon of Rayleigh Flat fading and its simulation using Clarke’s model and Young’s model were discussed in the previous posts. The performance (Eb/N0 Vs BER) of BPSK modulation (with coherent detection) over Rayleigh Fading channel and its comparison over AWGN channel is discussed in this post.

We first investigate the non-coherent detection of BPSK over Rayleigh Fading channel and then we move on to the coherent detection. For both the cases, we consider a simple flat fading Rayleigh channel (modeled as a – single tap filter – with complex impulse response – h). The channel also adds AWGN noise to the signal samples after it suffers from Rayleigh Fading.

The received signal y can be represented as

$$y=hx+n$$

where n is the noise contributed by AWGN which is Gaussian distributed with zero mean and unit variance and h is the Rayleigh Fading response with zero mean and unit variance. (For a simple AWGN channel without Rayleigh Fading the received signal is represented as y=x+n).

### Non-Coherent Detection:

In non-coherent detection, prior knowledge of the channel impulse response (“h” in this case) is not known at the receiver. Consider the BPSK signaling scheme with ‘x=+/- a’ being transmitted over such a channel as described above. This signaling scheme fails completely (in non coherent detection scheme), even in the absence of noise, since the phase of the received signal y is uniformly distributed between 0 and 2pi regardless of whether x[m]=+a or x[m]=-a is transmitted. So the non coherent detection of the BPSK signaling is not a suitable method of detection especially in a Fading environment.

### Coherent Detection:

In coherent detection, the receiver has sufficient knowledge about the channel impulse response.Techniques like pilot transmissions are used to estimate the channel impulse response at the receiver, before the actual data transmission could begin. Lets consider that the channel impulse response estimate at receiver is known and is perfect & accurate.The transmitted symbols (‘x’) can be obtained from the received signal (‘y’) by the process of equalization as given below.

$$\hat{y}=\frac{y}{h}=\frac{hx+n}{h}=x+z$$

here z is still an AWGN noise except for the scaling factor 1/h. Now the detection of x can be performed in a manner similar to the detection in AWGN channels.

The input binary bits to the BPSK modulation system are detected as

$$\begin{matrix} r=real(\hat{y})=real(x+z) \\ \; \; \hat{d} =1, \; \; if \;r> 0 \\ \; \; \hat{d}=0 , \; \; if \; r< 0 \end{matrix}$$

### Theoretical BER:

The theoretical BER for BPSK modulation scheme over Rayleigh fading channel (with AWGN noise) is given by

$$P_{b} =\frac{1}{2} \left ( 1-\sqrt{\frac{E_{b}/N_{0}}{1+E_{b}/N_{0}}}\right)$$

The theoretical BER for BPSK modulation scheme over an AWGN channel is given here for comparison

$$P_{b}=\frac{1}{2}erfc(\sqrt{E_{b}/N_{0}})$$

### Simulation Model:

The following model is used for the simulation of BPSK over Rayleigh Fading channel and its comparison with AWGN channel

### Matlab Code:

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

### Simulation Results:

The Simulated and theoretical performance curves (Eb/N0 Vs BER) for BPSK modulation over Rayleigh Fading channel and the AWGN is given below.

# BER Vs Eb/N0 for QPSK modulation over AWGN

(2 votes, average: 5.00 out of 5)

This post is a part of the ebook : Simulation of digital communication systems using Matlab – available in both PDF and EPUB format.

In the previous article we saw about how QPSK modulation and demodulation can be done. This concept is extended further to simulate the performance of QPSK modulation technique over an AWGN.

### Transmitter:

For the QPSK modulation , a series of binary input message bits are generated. In QPSK, a symbol contains 2 bits. The generated binary bits are combined in terms of two bits and QPSK symbols are generated. From the constellation of QPSK modulation the symbol ’00’ is represented by 1, ’01’ by j (90 degrees phase rotation), ’10’ by -1 (180 degrees phase rotation) and ’11’ by -j (270 degrees phase rotation). In pi/4 QPSK, these phase rotations are offset by 45 degrees. So the effective representation of symbols in pi/4-QPSK is ’00’=1+j (45 degrees), ’01’=-1+j (135 degrees), ’10’ = -1-j (225 degrees) and ’11’= 1-j (315 degrees).

Here we are simulating a pi/4 QPSK system.Once the symbols are mapped, the power of the QPSK modulated signal need to be normalized by $$\frac{1}{\sqrt{2}}$$.

### AWGN channel:

For QPSK modulation the channel can be modeled as

$$y=ax+n$$

where y is the received signal at the input of the QPSK receiver, x is the complex modulated signal transmitted through the channel , a is a channel amplitude scaling factor for the transmitted signal usually 1. ‘n’ is the Additive Gaussian White Noise random random variable with zero mean and variance $$\sigma^2$$. For AWGN the noise variance in terms of noise power spectral density $$N_0$$ is given by,

$$\sigma^{2}=\frac{N_{0}}{2}$$

For M-PSK modulation schemes including BPSK, the symbol energy is given by

$$E_s = R_m R_c E_b$$

where $$E_s$$ =Symbol energy per modulated bit ($$x$$), $$Rm = log_2(M)$$ , (for BPSK M=2, QPSK M=4, 16 QAM M=16 etc..,). $$R_c$$ is the code rate of the system if a coding scheme is used. In our case since no coding scheme is used Rc = 1. Eb is the Energy per information bit.

Assuming $$E_s=1$$ for BPSK (Symbol energy normalized to 1) $$\frac{E_b}{N_0}$$ can be represented as (using above equations),

$$\frac{E_b}{N_0} = \frac{E_s}{R_m R_c N_0}$$

$$\frac{E_b}{N_0} = \frac{E_s}{R_m R_c N_0} = \frac{E_s}{R_m R_c 2 \sigma^{2}} = \frac{1}{2 R_m R_c \sigma^{2}}$$

From the above equation the noise variance for the given $$\frac{E_b}{N_0}$$ can be calculated as

$$\sigma^{2} = \left( 2 R_m R_c \frac{E_b}{N_0}\right) ^{-1}$$

For the channel model randn function in Matlab is used to generate the noise term. This function generates noise with unit variance and zero mean. In order to generate a noise with sigma $$\sigma$$ for the given $$\frac{E_b}{N_0}$$ ratio , use the above equation , find $$\sigma$$, multiply the ‘randn’ generated noise with this sigma , add this final noise term with the transmitted signal to get the received signal. For a pi/4 QPSK system, since the modulated signal is in complex form, the noise should also be also in complex form.

QPSK receiver employs two threshold detectors that detect real(inphase arm) and imaginary part (quadrature arm). The detected signals are sent through a parallel to serial converter (achieved by “reshape” function in MATLAB).

The Eb/N0 Vs BER Curve for QPSK and BPSK are identical. This is because, QPSK actually consists of two orthogonal BPSK systems. Since the two individual BPSK systems are orthogonal, they don’t interfere with each other. That is why the BER curved for QPSK and BPSK are identical.

### Matlab Code:

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

# BER Vs Eb/N0 for BPSK modulation over AWGN

(7 votes, average: 3.57 out of 5)

In the previous article we saw about how Passband BPSK modulation and demodulation can be done. This concept is extended further to simulate the performance of BPSK modulation technique over an AWGN. Note that this is a baseband simulation of BPSK modulation and demodulation. Baseband simulation are faster and yields performace results same as that of pass band simulation.

### Transmitter:

For the BPSK modulation , a series of binary input message bits are generated of which ‘1’s are represented by 1v and ‘0’s are translated as ‘-1’ v (equivalent to NRZ coding as discussed in the previous post).

### AWGN channel:

For BPSK modulation the channel can be modeled as

$$y=ax+n$$

where y is the received signal at the input of the BPSK receiver, x is the modulated signal transmitted through the channel , a is a channel amplitude scaling factor for the transmitted signal usually 1. ‘n’ is the Additive Gaussian White Noise random random variable with zero mean and variance σ2.

For AWGN the noise variance in terms of noise power spectral density (N0) is given by,

$$\sigma^{2}= \frac{N_0}{2}$$

For M-ARY modulation schemes like M-PSK including BPSK, the symbol energy is given by,

$$E_s = R_m R_c E_b$$

where Es =Symbol energy per modulated bit (x), Rm = log2(M) , (for BPSK M=2, QPSK M=4, 16 QAM M=16 etc..,). Rc is the code rate of the system if a coding scheme is used. In our case since no coding scheme is used Rc = 1. Eb is the Energy per information bit.

Assuming Es=1 for BPSK (Symbol energy normalized to 1) Eb/N0 can be represented as (using above equations),

$$\frac{E_b}{N_0}=\frac{E_s}{R_m R_c N_0}$$

$$\frac{E_b}{N_0}=\frac{E_s}{R_m R_c N_0}=\frac{E_s}{R_m R_c 2 \sigma^{2}} = \frac{1}{2 R_m R_c \sigma^{2}}$$

From the above equation the noise variance for the given Eb/N0 can be calculated as

$$\sigma^{2} = \left( 2 R_m R_c \frac{E_b}{N_0}\right)^{-1}$$

For the channel model randn function in Matlab is used to generate the noise term. This function generates noise with unit variance and zero mean. In order to generate a noise with sigma σ for the given Eb/N0 ratio , use the above equation , find σ, multiply the ‘randn’ generated noise with this sigma , add this final noise term with the transmitted signal to get the received signal.

BPSK receiver can be a simple threshold detector which categorizes the received signal as ‘0’ or ‘1’ depending on the threshold that is being set. Calculation of Theoretical BER for BPSK over AWGN is discussed here.

### Matlab code:

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

### Simulation Result:

(7 votes, average: 3.57 out of 5)

# BPSK modulation and Demodulation

(4 votes, average: 4.50 out of 5)

This post is a part of the ebook : Simulation of digital communication systems using Matlab – available in both PDF and EPUB format.

### BPSK Modulation:

In digital modulation techniques a set of basis functions are chosen for a particular modulation scheme.Generally the basis functions are orthogonal to each other. Basis functions can be derived using ‘Gram Schmidt orthogonalization‘ procedure.Once the basis function are chosen, any vector in the signal space can be represented as a linear combination of the basis functions.

In Binary Phase Shift Keying (BPSK) only one sinusoid is taken as basis function modulation. Modulation is achieved by varying the phase of the basis function depending on the message bits. The following equation outlines BPSK modulation technique.

$$\begin{array}{l} S_{0}(t) = A cos (\omega t) \rightarrow \mbox{represents ‘0’}\\ S_{1}(t) = A cos (\omega t+\pi ) \rightarrow \mbox{represents ‘1’}\\ \end{array}$$

The constellation diagram of BPSK will show the constellation points lying entirely on the x axis. It has no projection on the y axis. This means that the BPSK modulated signal will have an in-phase component (I) but no quadrature component (Q). This is because it has only one basis function.

A BPSK modulator can be implemented by NRZ coding the message bits (1 represented by +ve voltage and 0 represented by -ve voltage) and multiplying the output by a reference oscillator running at carrier frequency ω.

### BPSK Demodulation:

For BPSK demodulator , a coherent demodulator is taken as an example. In coherent detection technique the knowledge of the carrier frequency and phase must be known to the receiver. This can be achieved by using a Costas loop or a PLL (phase lock loop) at the receiver. A PLL essentially locks to the incoming carrier frequency and tracks the variations in frequency and phase. For the following simulation , neither a PLL nor a Costas loop is used but instead we simple use the output of the PLL or Costas loop. For demonstration purposes we simply assume that the carrier phase recovery is done and simply use the generated reference frequency at the receiver (cos(ωt)).

In the demodulator the received signal is multiplied by a reference frequency generator (assuming the PLL/Costas loop be present). The multiplied output is integrated over one bit period using an integrator. A threshold detector makes a decision on each integrated bit based on a threshold. Since an NRZ signaling format is used with equal amplitudes in positive and negative direction, the threshold for this case would be ‘0’.

### File 1: NRZ_Encoder.m

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

### File 2: BPSK_Mod_Demod.m

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

### Simulation Results

The output of the Matlab code provides more insight into the modulation techniques. Apart from plotting the modulated and demodulated signal it also shows the constellation at transmitter/receiver and Power Spectral Density of the BPSK modulated signal.

# Colored Noise Generation in Matlab

(2 votes, average: 1.50 out of 5)

In communication systems the noise is mostly modeled as white noise. When there exists a noise that is “white” , then there must also exist a noise that is colored too.

White noise has constant power spectral density across the entire frequency spectrum (extending upto infinity).There is no correlation between the samples of a white noise process at different time instances i.e. the auto correlation or the auto covariance of white noise is zero for all lags except for lag L=0. So the auto covariance of white noise process will be an impulse function at lag L=0.

### Colored Noise:

When the power spectral density of the noise is not uniform across the entire frequency spectrum, it is called colored noise. For a comprehensive definition of various types of colored noise refer [1]. There exist non zero values for auto correlation or auto covariance at different time instances for the colored noise. The auto covariance is maximum for zero lag (L=0) and decreases gradually for increasing and decreasing values of lag (L).

### Generation of Colored Noise:

Colored noise can be generated by passing the white noise through a shaping filter. The shaping filter is a dynamic filter , usually a low pass filter. The response of the colored noise can be varied by adjusting the parameters of the shaping filter. The low pass filter can be implemented in various ways in Matlab. The Matlab’s “FILTER” function is used in this simulation.

Often it is stated in channel detector applications that a matched filter is the optimum detector “in the presence of white noise ” since it can maximize the output SNR. In the following posts simulation of optimum matched filter in the presence of white noise/colored noise will be demonstrated.

### Matlab code:

Check this book for full Matlab code.
Simulation of Digital Communication Systems Using Matlab – by Mathuranathan Viswanathan

### Sample Output:

(2 votes, average: 1.50 out of 5)