Download All

No BSD License  

Highlights from
Developing custom modulation schemes

image thumbnail
from Developing custom modulation schemes by Marta
Performing Custom Modulation and Demodulation using Agilent Instruments, MATLAB, and Simulink

Developing custom modulation schemes

Developing custom modulation schemes

This script shows how MATLAB can be used, together with Agilent instruments, for prototyping custom modulation schemes. It demonstrates MATLAB functionality for data acquisition, visualisation and analysis. It finishes by a simple example illustrating the first step towards simulate the prototype communications system in Simulink.

More precisely, this script shows how to:

  • Use MATLAB to prototype a communication system that consists of a random bit generator, bit-to-integer converter, baseband modulator with custom modulation format, channel, demodulator, integer-to-bit converter.
  • Compute the system's bit error rate (BER)
  • Upload generated custom waveform to an Agilent ESG Vector Signal Generator
  • Acquire the signal using Agilent MXA Signal Analyser
  • Demodulate and analyse the acquired signal.
  • Design a prototype Simulink model simulating a communication system using custom modulation scheme

The demodulation of the acquired signal is based on IS136 demo written by Todd Atkins and the functions for creating Agilent logo constellation and uploading the data on Agilent vector signal generator are from Tom Gaudette.

Contents

I. Prototype the communication system

Use MATLAB and Communications Toolbox to prototype the communications system using custom modulation scheme.

Set system parameters

Set the values of the main parameters

% size of signal constellation:     32
M = 32;
% number of bits per symbol:        5
k = log2(M);
% number of bits to process:        30000
n = 3e4;
% oversampling factor:              8
nSamp = 8;

Create signal source

Create a binary data stream as a column vector. Then plot first 40 bits in a stem plot.

% random binary data stream
x = randint(n,1);

stem(x(1:40), 'filled');
title('Random bits', 'FontSize', 18);
xlabel('Bit index', 'FontSize', 18);
ylabel('Binary value', 'FontSize', 18);

Map bits to symbols

Convert the bits in x into k-bit symbols. 7500 symbols (30000/4) ranging from 0 to 31

xsym = bi2de(reshape(x, k, length(x)/k).', 'left-msb');

% plot first 10 symbols in a stem plot.
figure;
stem(xsym(1:10));
title('Random symbols', 'FontSize', 18);
xlabel('Symbol index', 'FontSize', 18);
ylabel('Integer value', 'FontSize', 18);

Create desired constellation

MATLAB provides a lot of flexibility in generating the modulation schemes. Let's create a custom constellation, for example in shape of Agilent logo. 

% constellation is represented by 32 complex numbers
constellation = agilentlogo();

f = figure;
plot(constellation, '.', 'MarkerSize', 16);

% or create more customised plot...
s = repmat([2^8 2^7 2^6  2^2], [1,8]);
scatter(imag(constellation), real(constellation), s, 'filled');
axis off
set(f, 'Color', 'w')

Modulate the signal

Modulate with 32 general QAM.

% complex signal, 6000 symbols.
ytx = genqammod(xsym, constellation);

Perform the pulse shaping

Modulation is often followed by pulse shaping, and demodulation is often preceded by a filtering or an integrate-and-dump operation. This section presents an example involving rectangular pulse shaping. Rectangular pulse shaping repeats each output from the modulator a fixed number of times to create an up-sampled signal. Rectangular pulse shaping can be a first step or an exploratory step in algorithm development, though it is less realistic than other kinds of pulse shaping (e.g. using Raised Cosine Filter). 

[ytxs] = rectpulse(ytx, nSamp);

Apply channel model

Send signal over an AWGN channel.

EbNo = 10; % In dB
snr = EbNo + 10*log10(k) - 10*log10(nSamp);  % 16.021 db
yrxs = awgn(ytxs, snr, 'measured');

Downsample the signal

If the transmitter up-samples the modulated signal, then the receiver should down-sample the received signal before demodulating. The "integrate and dump" operation is one way to down-sample the received signal. 

yrx = intdump(yrxs, nSamp);

View constellation with a scatter plot

Create scatter plot of noisy signal and transmitted signal on the same axes. plotting 5000 symbols

h = scatterplot(yrx(1:5e3), 1, 0, 'b.');
hold on;
scatterplot(ytx(1:5e3), 1, 0, 'r.', h);
title('Received Signal');
legend('Received signal', 'Signal constellation');
% set axis ranges.
axis([-1.2, 1.2, -1.2, 1.2]);
hold off;

Demodulate the signal

Demodulate the signal using the same Modem object as for the modulation

zsym = genqamdemod(yrx, constellation);

Map symbols to bits

Undo the bit-to-symbol mapping performed earlier.

% convert integers to bits.
z = de2bi(zsym, 'left-msb');
% convert z from a matrix to a vector.
z = reshape(z.', numel(z), 1);

Computate BER

Compare x and z to obtain the number of errors and the bit error rate.

[number_of_errors, bit_error_rate] = biterr(x, z)

number_of_errors =

   801


bit_error_rate =

    0.0267

II. Test the system on real data using instruments

We have seen above how to use MATLAB and Communications Toolbox to prototype and test a communication system using custom modulation. In order to test the behaviour of the system on real data, we can replace the simulated channel by instruments, such as arbitrary waveform generator and signal analyser. In order to control the instruments and acquire the data we are going to use the Instrument Control Toolbox.

Upload data to ESG Vector Signal Generator

Upload the simulated signal to a Vector Signal Generator. In this example we have used AGILENT E4438C Vector Signal Generator. We uploaded the function using Agilent Waveform Download Assistant. 

% define the instrument host name
esgIp = 'E4438C_A16616';

% define the instrument sample rate
esgSampleRate = 22.5e6;

% upload the waveform and play it with sampleRate 22.5 MHz
% ESG_DL(ytxs, esgSampleClock, esgIp, 'SOURce:FREQuency 2E9', 'POWer -20')

% Here is the description of function ESG_DL
help('ESG_DL')

  ESG_DL allows the user to pass ESG control parameters via Simulink 'ESG Download Block'
  and then download the simulated data to the Agilent E4438C Signal Generator
  Function variables which get passed from Simulink are:
  'data'                Complex data structure generated from Simulink
  'Fs'                  Sample Rate
  'ipaddress'           IP Address of ESG
  'Source_Freq'         ESG RF Frequency
  'RF_Amplitude_Out'    ESG RF Amplitude 
  NOTE: The called functions below require the Agilent Waveform Download Assistant .m files loaded in your path 
  Written by David L. Barner - 12/7/2006
  Revision Marta Wilczkowiak - 01/10/2007

Receive data from MXA Signal Analyser

Acquire the IQ data. In this example we have used Agilent MXA Signal Analyser and MATLAB Instrument Driver for this instrument from MATLAB Central. To acquire the signal to MATLAB, you can use TMTOOL or one of the Agilent MXA examples from www.agilent.com/find/sa_programming. In this example the data was acquired from MXA Signal Analyzer in IQ analyzer mode with bandwith 25MHz and time span 2ms. The resulting MXA sample time was 22.221e-9s. This means that for each uploaded sample we acquired around 2 samples. 

% define the instrument IP address
mxaIp = '172.16.27.192';

% Acquire IQ vector
% iq = acquireIqVector(mxaIp);
% Here is the description of function acquireIqVector
help('acquireIqVector')

% compute the rate ration between the instruments
mxaSampleTime = 22.221e-9;
mxaSampleRate = 1/mxaSampleTime;
nInstrSamp = mxaSampleRate/esgSampleRate

  MATLAB/MXA example 7
  Getting IQ data using the MXA driver and plot display


nInstrSamp =

    2.0001

Save the IQ vector

Save the IQ vector to file acquiredIQ.mat.

% save acquiredIQ iq

Analyse the acquired signal

Use MATLAB algorithm development environment to locate the symbols in the acquired IQ vector.

Load IQ data

If you do not have hardware, you can perform the following analysis on the data saved in MAT file.

load acquiredIQ;

figure
plot(iq, 'r.')
title('Acquired IQ vector')
xlabel('In-Phase')
ylabel('Quadrature')

Rescale the data

Rescale the acquired data so that IQ values are within [-1, 1] range. The scaling is required only for the optional automatic step at the end. 

iq = (iq.');
maxVal = max( [ real( iq(:) ); imag( iq(:) ) ] );
iq = iq /(maxVal);

Isolate the Symbols

Because of differences between the rate at which symbols are broadcast, the rate at which samples are taken, and the blurring effect of transmission, the symbols are uniformly spaced throughout the acquired IQ data. The remainder of the samples are transitions between symbols. We know the symbol spacing from how signal is broadcasted and how it was acquired. By using the value manipulation tools in the editor and Cell mode, we can immediately visualize the results of different symbol choices.

In this step we allow also manual scaling of the signal, so that the symbol clusters most distant from the center lie on the unit circle. 

% manipulate the values below in Cell mode in order to separate the symbols
scale = 1.2;
s0 = 2;
ds = round(nSamp*nInstrSamp);

% subsample and scale the signal.
plot(iq(s0:ds:end)*scale, 'r.')
title('Subsampled IQ vector')
xlabel('In-Phase')
ylabel('Quadrature')

Extract symbols

Extract the symbols for the further processing

manualSymbols= iq(s0:ds:end)*scale;

Perform phase correction

At this point, we have located our symbols in the IQ data. However, the symbols still don't cluster into eight nice neat groupings. The problem is that we need to account for skew in the signal. If everything is not synchronized correctly, the symbols appear to "walk" around between their intended locations. We can correct for that with a linear phase correction. 

% modify this with cell mode (by .1 or less).
phaseValue = -0.15;

% create vector of complex numbers representing the increasing rotation
skewCorrection = exp(2*pi*i*linspace(0,phaseValue,length(manualSymbols)));
% multiply the original signal
manualSymbolsDeskewed = manualSymbols.*skewCorrection;

plot(manualSymbolsDeskewed, 'r.');
title('Subsampled IQ vector after phase correction')
xlabel('In-Phase')
ylabel('Quadrature')

Perform global rotation

Now that we have located our symbols and corrected for skew between the clocks, if our clusters aren't aligned correctly globally (all shifted clockwise or counter clockwise) we can apply a global rotation. 

% modify this with cell mode (by .1 or less).
rotation = -0.025;

% add global rotation
symbols = manualSymbolsDeskewed.*exp(2*pi*i*rotation);
plot(symbols,'r*');
title('Subsampled IQ vector after rotation correction')
xlabel('In-Phase')
ylabel('Quadrature')

III. Use pattern search to separate the symbols

At this point we have obtained the symbol separation using Cell mode, a part of MATLAB algorithm development environment. There are ways of finding the values allowing for the symbol separation automatically, using optimization methods. The cost function, corresponding to the distance of the separated symbols from the ideal points on the constellation contains a lot of local minima and in order to use classic optimisation methods we would need to know the approximate solutions. We can try to use global optimization methods, such as pattern search from Genetic Algorithm and Direct Search Toolbox in order to find the parameter values automatically without need of accurate starting point.

Initiate pattern search

Define starting point and explore the corresponding solution

% define the starting point
% [firstSymbol, scaling, phase, rotation]
x0 = [2, 1, 0, 0];

% sample rate should be known from transmition/acquisition sample rates so
% no need to optimize
sampleRate = round(nSamp*nInstrSamp);

% different parameters to be optimised have quite different values,
% which can influence the optimisation. Let's scale them in order to work
% with parameters with similar orders of magnitude
s = [1 1 100 100];

% show the current solution
showSolutionDemod(x0, sampleRate, iq);

% show the initial value of the cost function
cost = costSymbolSeparation(x0.*s, sampleRate, iq)

cost =

  334.4100

Run pattern search

Use fundtion PATTERNSEARCH to separate symbols

% define paramterews for pattern search
os = psoptimset('CompletePoll', 'on', 'CompleteSearch', 'on', 'SearchMethod', 'MADSPositiveBasis2N');

% set parameters bounds
lb = [1 0.5 -2*pi*100 -pi/2*100];
ub = [15,1.5,2*pi*100,pi/2*100];
% perform search
[xn, fval] = patternsearch(@(x)costSymbolSeparation(x, sampleRate, iq), ...
    x0.*s, [], [], [], [], lb, ub, [], os);

% de-scale the solution
x = xn./s

% show the current value of the cost function
fval

figure
showSolutionDemod(x, sampleRate, iq);

Optimization terminated: mesh size less than options.TolMesh.

x =

    4.2929    1.2679   -0.1541   -1.1460


fval =

   21.1228

IV. Simulate the communication system

We have seen above how to use MATLAB functionality to prototype and explore the properties of a communications system. We have seen also how to use instruments to test this system on real data. As a next step, we could simulate the behaviour of this system over time using Simulink. Here is a simple model simulating the system developed on the beginning of this script. It generates random bit sequences, transmits them and verifies the bit error rate. 

% todo colors, channel
custommodulationagilentlogoSL

Elaborate the model

The model above is fairly simple and does not take into account issues such as symbol timing recovery or phase recovery, which are necessary for processing life data. If it was to be implemented as a part of a real system, consideration would also have to be given to optimisation of the code for the performance, solving fixed points settings and partitioning of the system into components. For this particular problem, we have not developed models taking taking into consideration these real-world implementation issues, but we could give examples on how to tackle these problems using existing GPS or FM models as examples. 

% see the logo for the last time...
f = figure;
constellation = agilentlogo();
s = repmat([2^8 2^7 2^6  2^2], [1,8]);
scatter(imag(constellation), real(constellation), s, 'filled');
axis off
set(f, 'Color', 'w')


Published with MATLAB® 7.5

Contact us at files@mathworks.com