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);
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);
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);
%% 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;
--- /dev/null
+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
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)
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
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;
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');
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;
-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;
projects / 4yp.git / commitdiff