From: Adrian Iain Lam Date: Thu, 8 Nov 2018 21:49:44 +0000 (+0000) Subject: Added adaptive CMA equalizer, "fixed" CD and phase noise X-Git-Url: https://adrianiainlam.tk/git/?p=4yp.git;a=commitdiff_plain;h=5e9be3c421c4d52da9df842548f421751fa294d4 Added adaptive CMA equalizer, "fixed" CD and phase noise Implemented CMA equalizer. It can be seen that it works well in equalizing the received signal when CD is small (i.e. when the CD compensation does not work well). Convergence is to be discussed. Previous version of CD fails miserably when CD is small but compensation is not used. This is due to the constellation being rotated by the CD simulation. Code was added to rotate it back. Previous version of phase noise has a lot of errors due to cycle slips when attempting to recover the phase using the Viterbi-Viterbi algorithm. Experimenting with different block sizes did not help much. A differential PSK was attempted which seems to work well, agreeing with supervisor's predictions on asymptotic behaviour, though incurring a penalty at lower SNRs. --- diff --git a/CD_AWGN.m b/CD_AWGN.m index 1ea6ad8..be7de8c 100644 --- a/CD_AWGN.m +++ b/CD_AWGN.m @@ -10,7 +10,7 @@ sps = 4; % samples per symbol fs = Rsym * sps; % sampling freq (Hz) Tsamp = 1 / fs; -t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; EbN0_db = 0:0.2:14; EbN0 = 10 .^ (EbN0_db ./ 10); @@ -25,20 +25,22 @@ EsN0_db = 10 .* log10(EsN0); plotlen = length(EbN0); ber = zeros(1, plotlen); +berNoComp = zeros(1, plotlen); +berAdapt = zeros(1, plotlen); +berMatlabAdapt = zeros(1, plotlen); data = randi([0 M - 1], numSymbs, 1); -modData = pskmod(data, M, 0, 'gray'); +modData = pskmod(data, M, pi / M, 'gray'); x = txFilter(modData, rolloff, span, sps); %% Simulate chromatic dispersion -D = 20; % ps / (nm km) +D = 17; % ps / (nm km) lambda = 1550; % nm z = 10; % km xCD = chromaticDispersion(x, D, lambda, z, Tsamp); xCD = normalizeEnergy(xCD, numSymbs, 1); - for i = 1:plotlen snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); noiseEnergy = 10 ^ (-snr / 10); @@ -48,13 +50,54 @@ for i = 1:plotlen yCDComp = CDCompensation(y, D, lambda, z, Tsamp); r = rxFilter(yCDComp, rolloff, span, sps); + rNoComp = rxFilter(y, rolloff, span, sps); %% normalize energy %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy); rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps); - demodData = pskdemod(rSampled, M, 0, 'gray'); + rNoCompSampled = rNoComp(sps*span/2+1:sps:(numSymbs+span/2)*sps); + + %% rotate rNoCompSampled to match original data + theta = angle(-sum(rNoCompSampled .^ M)) / M; + %% if theta approx +pi/M, wrap to -pi/M + if abs(theta - pi / M) / (pi / M) < 0.1 + theta = -pi / M; + end + theta + rNoCompSampled = rNoCompSampled .* exp(-j * theta); + + %% adaptive filter + adaptFilterOut = adaptiveCMA(rSampled); + + demodData = pskdemod(rSampled, M, pi / M, 'gray'); + demodNoComp = pskdemod(rNoCompSampled, M, pi / M, 'gray'); + demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray'); + %%demodMatlabAdapt = pskdemod(matlabEq, M, pi / M, 'gray'); [bitErrors, ber(i)] = biterr(data, demodData); + [bitErrors, berNoComp(i)] = biterr(data, demodNoComp); + [~, berAdapt(i)] = biterr(data, demodAdapt); + %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt); + + + if EbN0_db(i) == 12 + figure(1); + scatterplot(rSampled); + title('Constellation after CD compensation'); + %%scatterplot(modData); + %%title('Original constellation'); + scatterplot(rNoCompSampled); + title('Constellation without CD compensation'); + scatterplot(adaptFilterOut); + title('Constellation with CD compensation and adaptive filter'); + %scatterplot(matlabEq); + %title('Matlab equalizer'); + ber(i) + %berNoComp(i) + berAdapt(i) + berMatlabAdapt(i) + end + end figure(1); @@ -63,9 +106,13 @@ clf; %% Plot simulated results semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); hold on; +%%semilogy(EbN0_db, berNoComp, 'g', 'LineWidth', 2); +semilogy(EbN0_db, berAdapt, 'm', 'LineWidth', 1.4); +%%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4); theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); -legend({'CD + AWGN + CD compensation', 'AWGN only'}, 'Location', 'southwest'); +legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ... + 'Theoretical AWGN'}, 'Location', 'southwest'); title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation')); grid on; diff --git a/adaptiveCMA.m b/adaptiveCMA.m new file mode 100644 index 0000000..084eab4 --- /dev/null +++ b/adaptiveCMA.m @@ -0,0 +1,57 @@ +function adaptFilterOut = adaptiveCMA(rSampled) + %% adaptive filter + %% CMA + taps = 19; % ODD taps + hxx = zeros(taps, 1); + %% hxx: real indices -K, ..., 0, ..., K. K = floor(taps/2) + %% MATLAB indices 1 1+K taps + %% initialize hxx, hxx[0] = 1, hxx[k] = hxx[-k] = 0 + hxx(ceil(taps/2)) = 1; + + mu = 1e-3; + numSymbs = length(rSampled); + + %% Check average energy of symbols + rSampledUnitEnergy = normalizeEnergy(rSampled, numSymbs, 1); + + adaptFilterOut = zeros(numSymbs, 1); + converged = 0; + convergeCount = 0; + + for it = 1:numSymbs + if it <= (taps - 1) / 2; + xp = [zeros((taps - 1) / 2 - it + 1, 1); rSampledUnitEnergy(1:it + (taps - 1) / 2)]; + elseif it + (taps - 1) / 2 > numSymbs + xp = [rSampledUnitEnergy(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)]; + else + xp = rSampledUnitEnergy(it - (taps - 1) / 2 : it + (taps - 1) / 2); + end + + xout = sum(hxx .* xp); + ex = 1 - abs(xout) ^ 2; + + if abs(ex) < 1e-3 + convergeCount = convergeCount + 1; + else + convergeCount = 0; + end + if ~converged && convergeCount >= 10 + converged = 1 + it + end + + adaptFilterOut(it) = xout; + + hxx = hxx + mu * ex * xout * conj(xp); + end + +%{ + %% try MATLAB builtin equalizer + alg = cma(mu); + eqObj = lineareq(taps, alg); + eqObj.Weights((taps + 1) / 2) = 1; + rPadded = [rSampledUnitEnergy; zeros((taps - 1) / 2, 1)]; + matlabEq = equalize(eqObj, rPadded); + matlabEq = matlabEq((taps + 1) / 2 : end); +%} +end diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m index db25965..7d9e980 100644 --- a/chromaticDispersion1Signal.m +++ b/chromaticDispersion1Signal.m @@ -1,63 +1,85 @@ M = 4; -numSymbs = 1000; +numSymbs = 100000; -%% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html Rsym = 2.5e10; % symbol rate (sym/sec) span = 6; % Tx/Rx filter span rolloff = 0.25; % Tx/Rx RRC rolloff sps = 4; % samples per symbol - fs = Rsym * sps; % sampling freq (Hz) Tsamp = 1 / fs; -t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; - +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; data = randi([0 M - 1], numSymbs, 1); -modData = pskmod(data, M, 0, 'gray'); +modData = pskmod(data, M, pi / M, 'gray'); x = txFilter(modData, rolloff, span, sps); %% Simulate chromatic dispersion -D = 20; % ps / (nm km) +D = 17; % ps / (nm km) lambda = 1550; % nm -z = 1000; % km +z = 10; % km [xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp); xCD = normalizeEnergy(xCD, numSymbs, 1); +EbN0_db = 8; +snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps); +noiseEnergy = 10 ^ (-snr / 10); +%%y = awgn(xCD, snr, 'measured'); y = xCD; - yCDComp = CDCompensation(y, D, lambda, z, Tsamp); +%%yCDComp = y; + +r = rxFilter(yCDComp, rolloff, span, sps); +rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps); + +%% if no CD comp, then rotate constellation. Use: +%{ +theta = angle(-sum(rSampled .^ M)) / M; +%% if theta approx +pi/M, wrap to -pi/M +if abs(theta - pi / M) / (pi / M) < 0.1 + theta = -pi / M; +end +rSampled = rSampled .* exp(-j * theta); +%} + +rAdaptEq = adaptiveCMA(rSampled); %% Compare original signal and compensated signal -figure(1); +figure(101); +clf; +tsym = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); subplot(211); -plot(real(x(1:300))); +plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)), 'b'); hold on -plot(real(yCDComp(1:300))); -hold off +plot(t(1:length(x)), real(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r'); +plot(tsym, real(rAdaptEq), 'xg'); +hold off; title('Real part'); -legend('original', 'dispersion compensated'); +legend('original', 'dispersion compensated', 'CMA equalized samples'); +axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]); subplot(212); -plot(imag(x(1:300))); -hold on -plot(imag(yCDComp(1:300))); -hold off +plot(t(1:length(x)), imag(normalizeEnergy(x, numSymbs*sps, 1)), 'b'); +hold on; +plot(t(1:length(x)), imag(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r'); +plot(tsym, imag(rAdaptEq), 'xg'); +hold off; title('Imag part'); - - -r = rxFilter(yCDComp, rolloff, span, sps); -r = normalizeEnergy(r, numSymbs, 1); % Add noise energy if needed - -rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps); +axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]); scatterplot(modData); title('Constellation of original modulation'); scatterplot(rSampled); -title('Constellation of sampled received waveform'); +title('Constellation of matched filter output'); +scatterplot(rAdaptEq); +title('Constellation of adaptive filter output'); + +demodData = pskdemod(rSampled, M, pi / M, 'gray'); +demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray'); -demodData = pskdemod(rSampled, M, 0, 'gray'); +[~, ber] = biterr(data, demodData) +[~, ber] = biterr(data, demodAdapt) diff --git a/phaseNoiseCorr.m b/phaseNoiseCorr.m index 1f4ed45..6d27f74 100644 --- a/phaseNoiseCorr.m +++ b/phaseNoiseCorr.m @@ -6,7 +6,10 @@ function [rPhaseEq, phiests] = phaseNoiseCorr(r, M, blocksize) block = r(l : min(l + blocksize - 1, length(r))); sum_M = sum(block .^ M); - phi_est = angle(sum_M) / M; % assume phase of 0 symbol is 0. + %% if phase of 0 symbol is 0, use: + phi_est = angle(sum_M) / M; + %% if phase of 0 symbol is pi/M, use: + %% phi_est = angle(-sum_M) / M; if l > 1 %% phase unwrapping diff --git a/phasenoise1signal.m b/phasenoise1signal.m index e0ee39d..fcd371a 100644 --- a/phasenoise1signal.m +++ b/phasenoise1signal.m @@ -2,14 +2,14 @@ numSymbs = 10000; M = 4; Rsym = 2.5e10; % symbol rate (sym/sec) -rolloff = 0.25; +rolloff = 0.5; span = 6; % filter span sps = 4; % samples per symbol fs = Rsym * sps; % sampling freq (Hz) Tsamp = 1 / fs; -t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; EbN0_db = 8; @@ -24,36 +24,33 @@ EsN0_db = 10 .* log10(EsN0); data = randi([0 M - 1], numSymbs, 1); -modData = pskmod(data, M, 0, 'gray'); +modData = dpskmod(data, M, 0, 'gray'); x = txFilter(modData, rolloff, span, sps); -linewidthTx = 0;%1e5; % Hz +linewidthTx = 0; % Hz linewidthLO = 1e6; % Hz -%%linewidthTx = Rsym * 1e-4; % Hz -%%linewidthLO = Rsym * 1e-3; % Hz [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp); snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps); -noiseEnergy = 10 ^ (-snr / 10); -y = awgn(xPN, snr, 'measured'); +%%y = awgn(xPN, snr, 'measured'); +y = xPN; r = rxFilter(y, rolloff, span, sps); -[rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps); -rPhaseEq = normalizeEnergy(rPhaseEq, numSymbs, 1 + noiseEnergy); -rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps); -demodData = pskdemod(rSampled, M, 0, 'gray')'; +rSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); +[rSaPhEq, phiests] = phaseNoiseCorr(rSampled, M, 40); -[bitErrors, ber] = biterr(data, demodData) +demodData = dpskdemod(rSaPhEq, M, 0, 'gray'); +[bitErrors, ber] = biterr(data, demodData.') figure(2); -plot(-phiests); +plot(repelem(-phiests, sps)); hold on; plot(pTxLO); legend('estimate', 'actual'); @@ -61,13 +58,14 @@ title('Phase noise estimation'); hold off; figure(3); -plot(t(1:length(x)), real(x)); +plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1))); hold on; -plot(t, real(rPhaseEq), 'r'); -%%sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); +%%plot(t, real(rPhaseEq), 'r'); +sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); %%plot(sampledTimes, real(rSampled), 'x'); +plot(sampledTimes, real(normalizeEnergy(rSaPhEq, numSymbs, 1)), 'x'); legend('original signal', 'corrected received signal'); -title('Phase noise correction, linewidth 1 MHz, E_b/N_0=8 dB'); +title('Phase noise correction, linewidth 1 MHz, $E_b/N_0=8$ dB'); ylabel('Real part of signals'); axis([t(1) t(300) -Inf +Inf]); hold off; diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m index a9c852a..d5185be 100644 --- a/phasenoise_AWGN.m +++ b/phasenoise_AWGN.m @@ -1,122 +1,101 @@ -function phasenoise_AWGN(rolloff, M, numSymbs) - %% Set defaults for inputs - if nargin < 3 - numSymbs = 1000; - end - if nargin < 2 - M = 2; - end - if nargin < 1 - rolloff = 0.25; - end - - plotted = 0; - - Rsym = 2.5e10; % symbol rate (sym/sec) - - span = 6; % filter span - sps = 4; % samples per symbol +numSymbs = 5e5; +M = 4; +rolloff = 0.5; - fs = Rsym * sps; % sampling freq (Hz) - Tsamp = 1 / fs; +Rsym = 2.5e10; % symbol rate (sym/sec) - t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; +span = 6; % filter span +sps = 4; % samples per symbol +fs = Rsym * sps; % sampling freq (Hz) +Tsamp = 1 / fs; - EbN0_db = 0:0.2:14; - EbN0 = 10 .^ (EbN0_db ./ 10); +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; - Es = 1; - Eb = Es / log2(M); - N0 = Eb ./ EbN0; - EsN0 = EbN0 .* log2(M); - EsN0_db = 10 .* log10(EsN0); +EbN0_db = 0:0.2:14; +EbN0 = 10 .^ (EbN0_db ./ 10); - plotlen = length(EbN0); +Es = 1; +Eb = Es / log2(M); +N0 = Eb ./ EbN0; - ber = zeros(1, plotlen); +EsN0 = EbN0 .* log2(M); +EsN0_db = 10 .* log10(EsN0); +plotlen = length(EbN0); - data = randi([0 M - 1], numSymbs, 1); - modData = pskmod(data, M, 0, 'gray'); +ber = zeros(1, plotlen); - x = txFilter(modData, rolloff, span, sps); +data = randi([0 M - 1], numSymbs, 1); +%%modData = dpskmod(data, M, pi / M, 'gray'); +modData = dpskmod(data, M, 0, 'gray'); - linewidthTx = 0;%1e5; % Hz - linewidthLO = 1e6; % Hz - %%linewidthTx = Rsym * 1e-4; % Hz - %%linewidthLO = Rsym * 1e-3; % Hz +x = txFilter(modData, rolloff, span, sps); - totalIterations = 1; +linewidthTx = 0; % Hz +%%linewidthLO = 1e6; % Hz +linewidthLO = Rsym * 1e-3; - for iter = 1:totalIterations - [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp); - for i = 1:plotlen - snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); - noiseEnergy = 10 ^ (-snr / 10); +avgSa = 40; - y = awgn(xPN, snr, 'measured'); - r = rxFilter(y, rolloff, span, sps); - %% normalize energy - %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy); +[xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp); - [rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps); +for i = 1:plotlen + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); + noiseEnergy = 10 ^ (-snr / 10); - rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps); - demodData = pskdemod(rSampled, M, 0, 'gray')'; + y = awgn(xPN, snr, 'measured'); - %%[bitErrors, ber(i)] = biterr(data, demodData); - [zzz, thisBER] = biterr(data, demodData); - ber(i) = ber(i) + thisBER / totalIterations; + r = rxFilter(y, rolloff, span, sps); + %% normalize energy + r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy); + rSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); + [rSaPhEq, phiests] = phaseNoiseCorr(rSampled, M, avgSa); - if plotted == 0 && EbN0_db(i) >6 && ber(i) > 1e-1 - plotted=1 - figure(1234); - plot(-phiests); - hold on; - plot(pTxLO); - legend('estimate', 'actual'); - hold off; + adaptiveFilterOut = adaptiveCMA(rSaPhEq.'); - figure(100); - %plot(t(1:length(x)), real(x)); - %%plot(t, real(x(1:length(t)))); - length(t) - length(x) - hold on - %%plot(t(1:length(xPhaseNoise)), real(xPhaseNoise)); + demodData = dpskdemod(rSaPhEq, M, 0, 'gray').'; + demodAdapt = dpskdemod(adaptiveFilterOut, M, 0, 'gray'); - plot(t, real(r), 'g') + [~, ber(i)] = biterr(data, demodData); - sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); + if EbN0_db(i) == 8 + figure(1234); + plot(repelem(-phiests, sps)); + hold on; + plot(pTxLO); + legend('estimate', 'actual'); + hold off; - plot(sampledTimes, real(rSampled), 'x') - hold off - - end - - end + figure(1); + scatterplot(rSaPhEq); + title('rSaPhEq'); end +end - figure(1); - clf; +figure(1); +clf; - %% Plot simulated results - semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); - hold on; +%% Plot simulated results +semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); +hold on; - theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); - legend({'Simulated phase noise', 'Without phase noise'}, 'Location', 'southwest'); +theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); +legend({'Simulated phase noise + correction', 'Theoretical AWGN'}, ... + 'Location', 'southwest'); - title(strcat(num2str(M), '-PSK with phase noise and correction')); - grid on; - xlabel('$E_b/N_0$ (dB)'); - ylabel('BER'); - formatFigure; -end +title({'QPSK with phase nosie and correction', ... + strcat(num2str(numSymbs * log2(M) / 1e3), '~kbit, LO~', ... + num2str(linewidthLO / 1e6), '~MHz, Av~', num2str(avgSa), ... + '~Sa')}); +grid on; +xlabel('$E_b/N_0$ (dB)'); +ylabel('BER'); + +formatFigure;