When a wireless signal travels from a transmitter (Tx) to a receiver (Rx) it undergoes some changes. In simple terms the signal s(t) is scaled by a factor h(t) and noise n(t) is added at the receiver. Let’s take this discussion forward with a simple example. Suppose the Tx transmits one of two possible symbols, +1 or -1. In technical lingo this is called Binary Phase Shift Keying (BPSK). If the channel scaling factor is 0.1 we will either get a +0.1 or -0.1 at the Rx to which AWGN noise is added. The noise is random in nature (having a Gaussian distribution) but for simplicity we assume that it can have one of two values, +0.01 or -0.01.
In a previous post we calculated the Bit Error Rate (BER) of a Massive MIMO system using two different channel models namely deterministic and probabilistic. The deterministic channel model is derived from the geometry of the array (ULA in this case) and the distribution of users in the cell. Whereas probabilistic channel model assumes that the channel is flat fading and can be modeled, between each transmit receive pair, as a complex, circularly symmetric, Gaussian random variable with mean of zero and variance of 0.5 per dimension.
In the previous two posts we discussed the fundamentals of array processing particularly the concept of beamforming (please check out array processing Part-1 and Part-2). Now we build upon these concepts to introduce some linear estimation techniques that are used in array processing. These are particularly suited to a situation where multiple users are spatially distributed in a cell and they need to be separated based upon their angles of arrival. But first let us introduce the linear model; I am sure you have seen this before.
x=Hs+w
Here, s is the vector of symbols transmitted by M users, H is the N x M channel matrix, w is the noise vector of length N and x is the observation vector of length N. The channel matrix formed by the channel coefficients is deterministic (as opposed to probabilistic) in nature as it is purely dependent upon the phase shifts that the channel introduces due to varying path lengths between the transmit and receive antennas. The impact of a channel coefficient can be thought of as a rotation of the complex signal without altering its amplitude.
This means that the channel acts like a single tap filter and the process of convolution is reduced to simple multiplication (a reasonable assumption if the symbol length is much larger than the channel delay spread). The channel model does not accommodate for path loss and fading that are also inherent characteristics of the channel. But the techniques are general enough for these effects to be factored in later. Furthermore, it is assumed that the channel H is known at the receiver. This is a realistic assumption if the channel is slowly varying and can be estimated by sending pilot signals.
Beamforming Using a Uniform Linear Array
So going back to the linear model we see that we know x and H while s and w are unknown. Here w cannot be estimated since it’s random in nature (remember what the term AWGN stands for?) and we ignore it for the moment. The structure of s is known. For example if we are using BPSK modulation then the m symbols of the signal vector s can either be +1 or -1. So we can start the process of symbol detection by substituting all possible combinations of s1, s2…sm and determine the combination that minimizes
||x-Hs||
This is called the Maximum Likelihood (ML) solution as it determines the combination that was most likely to have been transmitted based upon the observation.
Although ML is conceptually very appealing and yields good results it becomes prohibitively complex as the constellation size or number of transmit antennas increases. For example for 2-Transmit case and BPSK modulation there are 2^(1 bit x 2 antennas)=2^2=4 combinations, which seems quite simplistic. But if 16-QAM modulation is used and there are 4-Transmit antennas the number of combinations increases to 2^(4 bits x 4 antennas)=2^16=65536. So we conclude that ML is not the solution we are looking for if computational complexity is an issue (which might become less of an issue as the processing power of devices increases).
Next we turn our attention to a technique popularly known as Zero Forcing or ZF (the origins of the name I still do not know). According to this technique the channel has a multiplicative effect on the signal. So to remove this effect we simply divide the signal by the channel or in the language of matrices we perform matrix inversion. Mathematically we have:
x=Hs+w
H-1x=H-1Hs+H-1w
H-1x=s+H-1w=s̄
So we see that we get back the signal s but we also get a noise component enhanced by inverse of the channel matrix. This is the well-known problem of ZF called Noise Enhancement. Then there are other problems such as non-existence of the inverse when the channel H is not a square matrix (which only happens when the number of transmit and receive antennas is the same). The inverse of H also cannot be calculated if H is not full rank or determinant of H is zero. So we now introduce another technique called Least Squares (LS). According to this the signal vector can be estimated as
s̄=(HHH)-1HHx
This is also sometimes referred to as the Minimum Variance Unbiased Estimator, as described by Steven M. Kay in his classical book on Estimation Theory [Fundamentals of Statistical Signal Processing Vol-1]. This can be easily implemented in MATLAB using Moore Penrose Pseudo Inverse or pinv(H). This is much more stable than going for the direct inversion methods.
We next plot the Bit Error Rate (BER) using the code below. The number of receive antennas is varied from two to ten while the number of transmit antennas is fixed at four. The transmit antennas are assumed to be positioned at 30, 40, 50 and 60 degrees from the axis of the receive array. The receive antennas are separated by λ/2 meters. The frequency of operation is 1GHz but it is quite irrelevant to the scenario considered as everything is measured in multiples of wavelengths. The Eb/No ratio (roughly the signal to noise ratio) is varied from 5dB to 20dB in steps of 5dB.
Bit Error Rate for Changing Rx Array Length
As expected the BER for the two methods, other than ML, is more or less the same and decreases rapidly once the number of receive antennas becomes greater than number of transmit antennas (or number of signals). The case where the number of receive antennas is less than number of signals (equal powered and with a small angular separation) is dealt with by Overloaded Array Processing (OLAP) techniques and have been discussed in detail by James Hicks [Doctoral Dissertation] a student of Dr. Reed at Virginia Tech.
Strangely enough it is seen that the overloaded case is not the worst part of the BER curve. The worst BER is observed when the number the number of transmit and receive antennas is the same (four in this case). In other words the BER gradually increases as the rank of the channel matrix increases and then decreases once it reaches its maximum value. This is quite interesting and obviously has to do with Noise Enhancement that we discussed earlier. This will be further investigated in future posts.
For further information on the above methods visit this interesting article.
Update-1
So we struggled for a while to find out why the BER is worst at full rank and thought that there is something wrong in our model but ultimately we found that this has to do with how the pseudoinverse works and the way the tolerance limit (tol in MATLAB) for the singular values is set. We have found quite interesting results while experimenting with various inversion methods and the results are pending publication. Will keep you updated about the progress.
Update-2
We experimented with the MATLAB function pinv by changing the tolerance parameter. Previously we had used the default tolerance that is built into the function pinv. The default tolerance (tol) is defined as:
tol = max (size (H)) * sigma_max (H) * eps
where sigma_max (H) is the maximal singular value of channel matrix H
and eps is the machine precision.
More precisely, eps is the relative spacing between any two adjacent numbers in the machine’s floating point system. This number is obviously system dependent. On machines that support IEEE floating point arithmetic, eps is approximately 2.2204e-16 for double precision and 1.1921e-07 for single precision.
So back to the subject we experimented with two values of tol; 1.0 and 0.1 while changing the signal to noise ratio. The number of transmit antennas (users) is fixed at 4 while number of receive antennas is varied from 2 to 8. For tol value of 1.0 it is seen that changing the value of EbNo does not change the results much up to 6 receive antennas but after that the BER results rapidly diverge. For tol value of 0.1 the results are quite unexpected. The BER drops with increasing number of antennas up to N=5 but then there is an unexpected increase in the BER for N=6. This needs to be further investigated.
BER for tolerance of 1.0
BER for tolerance of 0.1
MATLAB CODE USED TO GENERATE ABOVE PLOT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MULTIUSER DETECTION USING
% A UNIFORM LINEAR ARRRAY
% AKA RECEIVE BEAMFORMING
% COPYRIGHT RAYMAPS (C) 2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
% SETTING THE PARAMETERS FOR THE SIMULATION
f=1e9; %Carrier frequency
c=3e8; %Speed of light
l=c/f; %Wavelength
d=l/2; %Rx array spacing
N=10; %Receive array length
theta=([30 40 50 60])*pi/180; %Angular placement of Tx array (users)
EbNo=10; %Energy per bit to noise PSD
sigma=1/sqrt(2*EbNo); %Standard deviation of noise
n=1:N; %Rx array vector
n=transpose(n); %Converting row to column
M=length(theta); %Tx array length
% RECEIVE SIGNAL MODEL (LINEAR)
s=2*(round(rand(M,1))-0.5); %BPSK signal of length M
H=exp(-i*(n-1)*2*pi*d*cos(theta)/l); %Channel matrix of size NxM
wn=sigma*(randn(N,1)+i*randn(N,1)); %AWGN noise of length N
x=H*s+wn; %Receive vector of length N
% PINV without tol
% y=pinv(H)*x;
% PINV with tol
y=pinv(H, 0.1)*x;
% DEMODULATION AND BER CALCULATION
s_est=sign(real(y)); %Demodulation
ber=sum(s!=s_est)/length(s); %BER calculation
As the data rates supported by wireless networks continue to rise the bandwidth requirements also continue to increase (although spectral efficiency has also improved). Remember GSM technology which supported 125 channels of 200KHz each, which was further divided among eight users using TDMA. Move on to LTE where the channel bandwidth could be as high as 20MHz (1.4MHz, 3MHz, 5MHz, 10MHz, 15MHz and 20MHz are standardized).
This advancement poses a unique challenge referred to as frequency selective fading. This means that different parts of the signal spectrum would see a different channel (different amplitude and different phase offset). Look at this in the time domain where the larger bandwidth means shorter symbol period causing intersymbol interference (as time delayed copies of the signal overlap on arrival at the receiver).
The solution to this problem is OFDM that divides the wideband signal into smaller components each having a bandwidth of a few KHz. Each of these components experiences a flat channel. To make the task of equalization simple a cyclic prefix (CP) is added in the time domain to make the effect of fading channel appear as circular convolution. Thus simplifying the frequency domain equalization to a simple division operation.
Shown below is the Python code that calculates the bit error rate (BER) of BPSK-OFDM which is the same as simple BPSK in a Rayleigh flat fading channel. However there is a caveat. We have inserted a CP which means we are transmitting more energy than simple BPSK. To be exact we are transmitting 1.25 (160/128) times more energy. This means that if this excess energy is accounted for the performance of BPSK-OFDM would be 1dB (10*log10(1.25)) worse than simple BPSK in Rayleigh flat fading channel.
Note:
Although we have shown the channel as a multiplicative effect in the figure above, this is only true for a single tap channel. For a multi-tap channel (such as the one used in the code above) the effect of the channel is that of a filter which performs convolution operation on the transmitted signal.
We have used a baseband model in our simulation and the accompanying figure. In reality the transmitted signal is upconverted before transmission by the antennas.
The above model can be easily modified for any modulation scheme such as QPSK or 16-QAM. The main difference would be that the signal would have a both a real part and an imaginary part, much of the simulation would remain the same. This would be the subject of a future post. For a MATLAB implementation of 64-QAM OFDM see the following post (64-QAM OFDM).
Serial to parallel and parallel to serial conversion shown in the above figure was not required as the simulation was done symbol by symbol (one OFDM symbol in the time domain represented 128 BPSK symbols in the frequency domain).
The channel model in the above simulation is quasi-static i.e. it remains constant for one OFDM symbol but then rapidly changes for the next, without any memory.
When wireless signals travel from a transmitter to a receiver they do so after reflection, refraction, diffraction and scattering from the environment. Very rarely is there a direct line of sight (LOS) between the transmitter and receiver. Thus multiple time delayed copies of the signal reach the receiver that combine constructively and destructively. In a sense the channel acts as an FIR (finite impulse response) filter. Furthermore since the transmitter or receiver may be in motion the amplitude and phase of these replicas varies with time.
There are several methods to model the amplitude and phase of each of these components. We look at one method called the “Smiths Fading Simulator” which is based on Clark and Gans model. The simulator can be constructed using the following steps.
1. Define N the number of Gaussian RVs to be generated, fm the Doppler frequency in Hz, fs the sampling frequency in Hz, df the frequency spacing which is calculated as df=(2*fm)/(N-1) and M total number of samples in frequency domain which is calculated as M=(fs/df).
2. Generate two sequences of N/2 complex Gaussian random variables. These correspond to the frequency bins up to fm. Take the complex conjugate of these sequences to generate the N/2 complex Gaussian random variables for the negative frequency bins up to -fm.
3. Multiply the above complex Gaussian sequences g1 and g2 with square root of the Doppler Spectrum S generated from -fm to fm. Calculate the spectrum at -fm and +fm by using linear extrapolation.
4. Extend the above generated spectra from -fs/2 to +fs/2 by stuffing zeros from -fs/2 to -fm and fm to fs/2. Take the IFFT of the resulting spectra X and Y resulting in time domain signals x and y.
5. Add the absolute values of the resulting signals x and y in quadrature. Take the absolute value of this complex signal. This is the desired Rayleigh distributed envelope with the required temporal correlation.
The Matlab code for generating Rayleigh random sequence with a Doppler frequency of fm Hz is given below.
Diffraction is a phenomenon where electromagnetic waves (such as light waves) bend around corners to reach places which are otherwise not reachable i.e. not in the line of sight. In technical jargon such regions are also called shadowed regions (the term again drawn from the physics of light). This phenomenon can be explained by Huygen’s principle which states that “as a plane wave propagates in a particular direction each new point along the wavefront is a source of secondary waves”. This can be understood by looking at the following figure. However one peculiarity of this principle is that it is unable to explain why the new point source transmits only in the forward direction.
Diffraction is Difficult to Model
The electromagnetic field in the shadowed region can be calculated by combining vectorially the contributions of all of these secondary sources, which is not an easy task. Secondly, the geometry is usually much more complicated than shown in the above figure. For example consider a telecom tower transmitting electromagnetic waves from a rooftop and a pedestrian using a mobile phone at street level. The EM waves usually reach the receiver at street level after more than one diffraction (not to mention multiple reflections). However, an approximation that works well in most cases is called knife edge diffraction, which assumes a single sharp edge (an edge with a thickness much smaller than the wavelength) separates the transmitter and receiver.
Knife Edge Model
The path loss due to diffraction in the knife edge model is controlled by the Fresnel Diffraction Parameter which measures how deep the receiver is within the shadowed region. A negative value for the parameter shows that the obstruction is below the line of sight and if the value is below -1 there is hardly any loss. A value of 0 (zero) means that the transmitter, receiver and tip of the obstruction are all in line and the Electric Field Strength is reduced by half or the power is reduced to one fourth of the value without the obstruction i.e. a loss of 6dB. As the value of the Fresnel Diffraction Parameter increases on the positive side the path loss rapidly increases reaching a value of 27 dB for a parameter value of 5. Sometimes the exact calculation is not needed and only an approximate calculation, as proposed by Lee in 1985, is sufficient.
Fresnel Diffraction Parameter (v) is defined as:
v=h√(2(d1+d2)/(λ d1 d2))
where
d1 is the distance between the transmitter and the obstruction along the line of sight
d2 is the distance between the receiver and the obstruction along the line of sight
h is the height of the obstruction above the line of sight
and λ is the wavelength
The electrical length of the path difference between a diffracted ray and a LOS ray is equal to φ=(π/2)(v²) and the normalized electric field produced at the receiver, relative to the LOS path is e-jφ. Performing a summation of all the exponentials above the obstruction (from v to positive infinity) gives us the Fresnel Integral, F(v).
Knife Edge Diffraction Model Using Huygens Principle
Plot of Diffraction Loss
The MATLAB codes used to generate the above plots are given below (approximate method followed by the exact method). Feel free to use them in your simulations and if you have a question drop us a comment.
MATLAB Code for Approximate Calculation of Diffraction Loss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculation of the path loss based on the value of % Fresnel Diffraction Parameter as proposed by Lee % Lee W C Y Mobile Communications Engineering 1985 % Copyright www.raymaps.com %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all v=-5:0.01:5; for n=1:length(v) if v(n) <= -1 G(n)=0; elseif v(n) <= 0 G(n)=20*log10(0.5-0.62*v(n)); elseif v(n) <= 1 G(n)=20*log10(0.5*exp(-0.95*v(n))); elseif v(n) <= 2.4 G(n)=20*log10(0.4-sqrt(0.1184-(0.38-0.1*v(n))^2)); else G(n)=20*log10(0.225/v(n)); end end plot(v, G, 'b') xlabel('Fresnel Diffraction Parameter') ylabel('Diffraction Loss (dB)')
MATLAB Code for Exact Calculation of Diffraction Loss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exact calculation of the path loss (in dB) % based on Fresnel Diffraction Parameter (v) % T S Rappaport Wireless Communications P&P % Copyright www.raymaps.com %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all v=-5:0.01:5; for n=1:length(v) v_vector=v(n):0.01:v(n)+100; F(n)=((1+1i)/2)*sum(exp((-1i*pi*(v_vector).^2)/2)); end F=abs(F)/(abs(F(1))); plot(v, 20*log10(F),'r') xlabel('Fresnel Diffraction Parameter') ylabel('Diffraction Loss (dB)')
We have used the following equations in the exact calculation of the Diffraction Loss [1] above. We did not want to scare you with the math so have saved it for the end.
Also please checkout this interesting video explaining the phenomenon of diffraction.
