Nerd Rage

A while back I posted an arti­cle on Sig­nal Sta­tis­tics for OFDM. Orig­i­nally, this would have included OTDM also, but it grew too long. So it was split into two parts. This is the sec­ond part. Please read the orig­i­nal arti­cle first, as it explains much of the ter­mi­nol­ogy and the present arti­cle by itself may be dif­fi­cult to fol­low.

OTDM Recap

Some basics of OTDM / Nyquist WDM are explained in The Art of Nyquist WDM which I rec­om­mend read­ing if you’re new to the con­cept. Recently, I find the name OTDM much more fit­ting due to its dual­ity with OFDM. How­ever, what’s com­monly termed OTDM in optics (mul­ti­plex­ing really short pulse trains) I assume stands for opti­cal TDM, since it doesn’t make any use of an orthog­o­nal­ity rela­tion (though tech­ni­cally, the trib­u­taries are orthog­o­nal).

Any­way, the basic premise of OTDM is that $K$ pulses are over­lap­ping in time. Only by their spe­cial shape (the trun­cated sinc func­tion) can they be detected with­out crosstalk. The sinc func­tion has the for­tu­nate prop­erty of hav­ing zeros at reg­u­lar inter­vals. By mak­ing it so that these zeros occur at the cen­ters of all neigh­bor­ing sym­bol slots, there will only be one sym­bol not equal to zero at each slot cen­ter (this being the sym­bol that is cen­tered on that par­tic­u­lar slot). At all other times within the sym­bol slot, all sym­bols can have val­ues other than zero and inter­fere. Hence, we imme­di­ately expect quite dif­fer­ent sig­nal sta­tis­tics depend­ing on the posi­tion within the sym­bol slot we are inter­ested in.

The spec­trum of the sinc func­tion is rec­tan­gu­lar (or almost rec­tan­gu­lar for the trun­cated sinc) so that many such chan­nels could be put very close to each other with­out spec­tral crosstalk.

Impor­tance Sam­pling

In con­trast to OFDM which was (in a spe­cial case) the super­po­si­tion of ran­domly cho­sen let­ters from the same alpha­bet (the QAM con­stel­la­tion dia­gram) with equal weights, the weights of the dif­fer­ent super­posed sym­bols in OTDM are defined by the sinc func­tion. At the sym­bol slot cen­ter, for instance, the weight for the sym­bol cen­tered in that slot is unity while for all other sym­bols it is zero. In this case the PMF for a quad­ra­ture, the sig­nal or its power can be sim­ply deter­mined from the con­stel­la­tion dia­gram of the mod­u­la­tion for­mat.

At all other posi­tions the weights are dif­fer­ent and non-zero. We can­not there­fore use com­bi­na­tor­i­cal means to tackle this prob­lem. Due to the low prob­a­bil­i­ties involved (which are expected to be sim­i­lar to OFDM), an exhaus­tive numer­i­cal sim­u­la­tion is not fea­si­ble and a Monte-Carlo sim­u­la­tion will not yield the rel­e­vant low-probability details. We can, how­ever, apply the con­cept of impor­tance sam­pling to our Monte Carlo sim­u­la­tions [1]. Here, we bias the ran­dom draw­ing (of points from the con­stel­la­tion dia­gram for each sym­bol) to favor the occur­rence of sig­nals with high pow­ers, i.e. we pre­fer to stay in the cor­ners of the con­stel­la­tion dia­gram and switch to the oppo­site cor­ner when­ever the sign of the par­tic­u­lar sinc-weight for a sym­bol becomes neg­a­tive. In this way we get sam­ples for which all the con­stituents are quite large and in phase. We have to retain some ran­dom­ness, how­ever, since we are inter­ested in the PMF/CDF and not just the max­i­mum pos­si­ble value (more on that below). The biased prob­a­bil­ity for the sym­bols is com­pen­sated when later cal­cu­lat­ing the PMF.

Bias­ing toward higher sam­ple pow­ers results in a smaller num­ber of sam­ples with low power, which in turn leads to larger vari­ance in that part of the PMF. We would like to have an accu­rate rep­re­sen­ta­tion of the whole PMF, though. Hence, we will gen­er­ally sim­u­late the same mod­u­la­tion for­mat with dif­fer­ent bias val­ues and later com­bine their results using the bal­ance heuris­tic which is basi­cally just a pretty opti­mal way of com­bin­ing such data [2].

You will find the MATLAB code (not opti­mized for speed or any­thing) at the end of the arti­cle.

Sig­nal Power

Since the code is avail­able to play around with, we will look at just a few results to high­light the dif­fer­ences between OFDM and OTDM. In par­tic­u­lar, the sym­bol length will be $16T$, i.e. 16 times the time between con­sec­u­tive sym­bols (or the inverse sym­bol rate), and thus $K=16$. The alpha­bets used will be 16-QAM ($N=4$) and 64-QAM ($N=8$). As can be seen in Fig. 5 in the OFDM arti­cle, this results in prob­a­bil­i­ties in the vicin­ity of $10^{-28}$ for the largest pow­ers in OFDM sig­nals, some­thing not fea­si­ble to be fig­ured out by “reg­u­lar” Monte Carlo tech­niques.

The fig­ure shows the worst case regard­ing high sig­nal pow­ers, with the sam­pling instants cho­sen at the sym­bol slot bound­aries. See also the next sec­tion for more details. While the prob­a­bil­i­ties of the max­i­mum val­ues are the same (they depend only on $K$ and $N$ as shown in the OFDM post), the max­i­mum observ­able power is less than half in OTDM. At high prob­a­bil­i­ties like $10^{-6}$ the dif­fer­ence is smaller. This is a result of the het­ero­ge­neous weights of the con­tri­bu­tions from the dif­fer­ent sym­bols to each sam­ple, as opposed to equal (mag­ni­tudes of) weigths in OFDM.

Some­times there will be “kinks” in the graphs so obtained by impor­tance sam­pling. These are a result of out­liers whose effect is ampli­fied by the algo­rithm (in the code, data with very low dataweights con­tribute sig­nif­i­cantly more). This is a known prob­lem of impor­tance sam­pling. As these occur ran­domly, they might or might not be there. Their posi­tion is also ran­dom, so that they can be iden­ti­fied by mul­ti­ple runs of the code.

Other Sam­ple Times

In OFDM, the sam­ple power dis­tri­b­u­tion var­ied only lit­tle with sam­ple tim­ing within the sym­bol slot (cf. Figs. 6 and 7 in the OFDM post). As the con­stel­la­tion dia­grams of con­sec­u­tive OTDM sym­bols do not rotate rel­a­tive to each other, it is quite straight­for­ward to deter­mine the max­i­mum pos­si­ble sam­ple power at var­i­ous instants within the sym­bol slot by sim­ply adding the weight mag­ni­tudes (plus squar­ing and scal­ing). They can even be given analytically.$^1$

Unlike OFDM, the max­i­mum power does not have mul­ti­ple local max­i­mums within the sym­bol slot, but varies by about an order of mag­ni­tude between sym­bol slot cen­ter and bound­aries. The sam­ple power CDFs vary accord­ingly:

For $t_S=0.5$ (the sym­bol slot cen­ter), the CDF cor­re­sponds to the con­stel­la­tion dia­gram, as at this instant only one of the over­lap­ping sym­bols is non-zero.

Over­sam­pling

For (always nec­es­sary) over­sam­pling, the dis­tri­b­u­tions for the appro­pri­ate sam­pling instances can sim­ply be com­bined. Inter­est­ingly, for sim­ple two-fold over­sam­pling, half the sam­ples will be taken in the sym­bol slot cen­ter and the other half at its bound­aries. Thus, the worst case regard­ing sym­bol power is included and the high-power sam­ples result­ing thereof need to be accounted for by leav­ing appro­pri­ate head­room in the digital-analog con­vert­ers, ampli­fiers and what­not.

MATLAB Code

clear variables

LoS = 16; % length of symbol in units of symbol slots
QAM = 64; % number of possible values in each quadrature for 16-QAM - only powers of 2 allowed
NoS = 2^20; % number of Monte Carlo samples

offset = 0.5; % relative sampling position within symbol slot - 0.5 is center

% list of biasing coefficients used; need not be integer
explist = [0,1,2,3]; % works well for 16-QAM, LoS = 16
% explist = [0,2,4,6,8]; % works well for 64-QAM, LoS = 16

bins = linspace(0,16,1001); % set histogram bins for graphing

%%% set-up of some variables

normalization = 1/sqrt(2/3*(QAM-1)); % alphabet normalization
normalization = normalization / 0.987345; % unit energy per symbol for LoS = 16
% normalization = normalization / 0.993669; % unit energy per symbol for LoS = 32
normweight = (1/QAM)^LoS; % probability of any particular symbol w/o biasing

% sinc coefficients of each sample within the symbol
coeffs = sinc((-LoS/2:LoS/2-1) + offset);

% constellation diagram
NoV = sqrt(QAM); % number of possible values in I and Q
values = 0:2:2*(NoV-1); values = values - NoV + 1; % possible values in I and Q for each symbol
values = values * normalization;
symbols = repmat(values, [NoV,1]) + 1i * repmat(values.', [1,NoV]);
symbols = symbols(:); % complex QAM alphabet

bins = [bins, Inf];
PMF = zeros(length(bins), length(explist));
bincounts = PMF;

%%% importance sampling Monte Carlo simulation

for jj = 1:length(explist)

    % biasing is based on distance to opposite corner symbols (depending on coeff sign)
    distplus = abs(symbols - symbols(1));
    wplus = 1 ./ (1 + distplus).^explist(jj); % weights avoid division by zero
    wplus = wplus / sum(wplus); % normalization

    distminus = abs(symbols - symbols(QAM)); % opposite corner
    wminus = 1 ./ (1 + distminus).^explist(jj);
    wminus = wminus / sum(wminus);

    origdata = rand(NoS,LoS); % generate random data (symbol plus neighbors)

    data = zeros(NoS,LoS);
    dataweights = zeros(NoS,LoS);

    for kk = 1:LoS
        if coeffs(kk) >= 0
            for ii = (QAM):-1:1; % draw from biased distribution
                data(origdata(:,kk) < sum(wplus(1:ii)),kk) = symbols(ii);
                dataweights(origdata(:,kk) < sum(wplus(1:ii)),kk) = wplus(ii);
            end
        else % for negative coeffs bias towards opposite corner
            for ii = (QAM):-1:1;
                data(origdata(:,kk) < sum(wminus(1:ii)),kk) = symbols(ii);
                dataweights(origdata(:,kk) < sum(wminus(1:ii)),kk) = wminus(ii);
            end
        end
    end

    data = data .* repmat(coeffs, [NoS, 1]); % each row corresponds to one sample

    powers = abs(sum(data, 2)).^2;
    dataweights = prod(dataweights, 2);

    if explist(jj) ~= 0 % remove IS outliers
        index = find(dataweights >= normweight);
        powers = powers(index);
        dataweights = dataweights(index) * NoS / 4;
            % division by 4 is necessary to compensate for biasing towards only 1 corner
    else
        dataweights = dataweights * NoS;
    end

    % build PMF histograms for each biasing value
    for ii = 1:length(bins)-1
        index = find(powers >= bins(ii) & powers < bins(ii+1));
        PMF(ii,jj) = sum(normweight ./ dataweights(index));
        bincounts(ii,jj) = length(index);
    end

end

% average histograms using balance heuristic
weights = PMF .* bincounts ./ repmat(sum(PMF .* bincounts, 2), [1, length(explist)]);
weights(isnan(weights)) = 0; % bincount of zero causes weights to become NaN

avgPMF = sum(PMF .* weights, 2);

% complementary CDF
for ii=1:length(avgPMF)
    CDF(ii) = sum(avgPMF(ii:end));
end

Note the nor­mal­iza­tion of the sym­bol ampli­tude in order to obtain unit mean energy per sym­bol. This accounts for the dis­tri­b­u­tion of the con­stel­la­tion points for the par­tic­u­lar QAM for­mat as well as the trun­cated sinc pulse shape. Since

$$\intop_{-K/2}^{K/2} \mathrm{sinc}^2 x \, dx < 1 \quad \text{for} \quad K < \infty$$

the energy reduc­tion result­ing from trun­ca­tion is com­pen­sated by appro­pri­ate scal­ing.

---

1 The main prob­lem is to find

$$S_\mathrm{max}^2 = \left[\sum_{k=0}^{K-1} \left| \mathrm{sinc} \left(k - K/2 + t_S\right) \right| \right]^2 = \left[\sum_{k=0}^{K-1} \left| \frac{\sin \bigl[\pi \left(k - K/2 + t_S\right)\bigr]}{\pi \left(k - K/2 + t_S\right)} \right| \right]^2$$

which is then mul­ti­plied by the high­est power occur­ing in the con­stel­la­tion dia­gram of the par­tic­u­lar mod­u­la­tion for­mat (which also depends on $K$ when prop­erly nor­mal­ized). For sim­plic­ity, we’ll assume $K$ to be even and we will only try and find $S_\mathrm{max}$, and leave the squar­ing to the inter­ested reader. The sine in the numer­a­tor will then have a con­stant value, deter­mined by $t_S$ that sim­ply alter­nates its sign with $k$ and thus can be writ­ten

$$S_\mathrm{max} = \frac{\sin \pi t_S}{\pi}\sum_{k=0}^{K-1} \left| \frac{c_S \left(-1\right)^k}{k - K/2 + t_S} \right|$$

with

$$c_S = \begin{cases} 1 & K/2 \text{ even}\\ -1 & K/2 \text{ odd}\end{cases}$$

We split the sum to get rid of the neg­a­tive denom­i­na­tor and with it the absolute value func­tion and the sign coef­fi­cients:

$$S_\mathrm{max} = \frac{\sin \pi t_S}{\pi}\sum_{k=0}^{K/2-1} \frac{1}{k + t_S} - \frac{1}{k - K/2 + t_S}$$

These are then just two gen­eral har­monic series. Their solu­tion can be given by gen­er­al­ized har­monic num­bers or in terms of the Digamma func­tion $\psi(\cdot)$ [3]

$$S_\mathrm{max} = \frac{\sin \pi t_S}{\pi} \left[ \psi\left(t_S - \frac{K}{2}\right) + \psi\left(t_S + \frac{K}{2}\right) - 2 \psi\left(t_S \right) \right]$$

[1] Impor­tance Sam­pling [Wikipedia]
[2] E. Veach, Robust Monte Carlo meth­ods for light trans­port sim­u­la­tion, dis­ser­ta­tion, Stan­ford Uni­ver­sity, 1998.
[3] T. M. Ras­sias and H. M. Sri­vas­tava, “Some classes of infi­nite series asso­ci­ated with the Rie­mann zeta and polygamma func­tions and gen­er­al­ized har­monic num­bers,” Applied Math­e­mat­ics and Com­pu­ta­tion, vol. 131, pp. 593–605, 2002.

---

last posts in Nyquist WDM:

• the Epic Fail of Nyquist WDM
• the Art of Nyquist WDM