Added presentation; DE-QPSK; CD with FFT; split-step Fourier

[4yp.git] / phasenoise_AWGN.m

diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m
index d5185be..f055bcc 100644 (file)
--- a/phasenoise_AWGN.m
+++ b/phasenoise_AWGN.m
@@ -1,4 +1,4 @@
-numSymbs = 5e5;
+numSymbs = 5e3;
 M = 4;
 rolloff = 0.5;
 
@@ -25,55 +25,91 @@ EsN0_db = 10 .* log10(EsN0);
 
 plotlen = length(EbN0);
 
-ber = zeros(1, plotlen);
+berPSK = zeros(1, plotlen);
+berDEPSK = zeros(1, plotlen);
+berDPSK = zeros(1, plotlen);
 
 data = randi([0 M - 1], numSymbs, 1);
-%%modData = dpskmod(data, M, pi / M, 'gray');
-modData = dpskmod(data, M, 0, 'gray');
 
-x = txFilter(modData, rolloff, span, sps);
-
-linewidthTx = 0; % Hz
-%%linewidthLO = 1e6; % Hz
-linewidthLO = Rsym * 1e-3;
-
-
-avgSa = 40;
-
-
-[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);
-
-  y = awgn(xPN, snr, 'measured');
-
-  r = rxFilter(y, rolloff, span, sps);
-  %% normalize energy
-  r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
+pskSym = pskmod(data, M, pi / M, 'gray');
+%% DEPSK: Part VII, M.G. Taylor (2009)
+depskSym = pskmod(data, M, 0, 'gray');
+for i = 2:numSymbs
+  depskSym(i) = depskSym(i) * depskSym(i-1);
+end
 
-  rSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-  [rSaPhEq, phiests] = phaseNoiseCorr(rSampled, M, avgSa);
+dpskSym = dpskmod(data, M, pi / M, 'gray');
 
-  adaptiveFilterOut = adaptiveCMA(rSaPhEq.');
+xPSK = txFilter(pskSym, rolloff, span, sps);
+xDEPSK = txFilter(depskSym, rolloff, span, sps);
+xDPSK = txFilter(dpskSym, rolloff, span, sps);
 
-  demodData = dpskdemod(rSaPhEq, M, 0, 'gray').';
-  demodAdapt = dpskdemod(adaptiveFilterOut, M, 0, 'gray');
+linewidthTx = 0; % Hz
+linewidthLO = 5e6; % Hz
+%linewidthLO = Rsym * 1e-3;
 
-  [~, ber(i)] = biterr(data, demodData);
+iterations = 25;
+avgSa = 40;
 
-  if EbN0_db(i) == 8
-    figure(1234);
-    plot(repelem(-phiests, sps));
-    hold on;
-    plot(pTxLO);
-    legend('estimate', 'actual');
-    hold off;
+for it = 1 : iterations
+  [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
+  [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
+  [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
+
+  for i = 1:plotlen
+    snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
+    noiseEnergy = 10 ^ (-snr / 10);
+
+    yPSK = awgn(xPSKpn, snr, 'measured');
+    yDEPSK = awgn(xDEPSKpn, snr, 'measured');
+    yDPSK = awgn(xDPSKpn, snr, 'measured');
+
+    rPSK = rxFilter(yPSK, rolloff, span, sps);
+    rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
+    rDPSK = rxFilter(yDPSK, rolloff, span, sps);
+
+    rPSKSamp = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+    rDEPSKSamp = rDEPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+    rDPSKSamp = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+
+    [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
+    [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
+
+    demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').';
+    %% The decoding method described in Taylor (2009)
+    %% works on the complex symbols, i.e. after taking
+    %% the nearest symbol in the constellation, but before
+    %% converting them back to integers/bits.
+    %% MATLAB's pskdemod() does not provide this intermediate
+    %% result, so to be lazy, a pskmod() call is performed
+    %% to obtain the complex symbols.
+    demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').';
+    remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray');
+    delayed = [1; remodDEPSK(1:end-1)];
+    demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray');
+
+    demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');
+
+    [~, ber] = biterr(data, demodPSK);
+    berPSK(i) = berPSK(i) + ber / iterations;
+    [~, ber] = biterr(data, demodDEPSK);
+    berDEPSK(i) = berDEPSK(i) + ber / iterations;
+    [~, ber] = biterr(data, demodDPSK);
+    berDPSK(i) = berDPSK(i) + ber / iterations;
+
+    if EbN0_db(i) == 8 && it == 1
+      figure(1234);
+      plot(repelem(-phiestsPSK, sps));
+      hold on;
+      plot(pTxLoPSK);
+      legend('estimate', 'actual');
+      hold off;
+
+      figure(1);
+      scatterplot(rPSKSampEq);
+      title('rPSKSampEq');
+    end
 
-    figure(1);
-    scatterplot(rSaPhEq);
-    title('rSaPhEq');
   end
 end
 
@@ -82,18 +118,30 @@ figure(1);
 clf;
 
 %% Plot simulated results
-semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
+semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
 hold on;
+semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
+semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);
 
 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-legend({'Simulated phase noise + correction', 'Theoretical AWGN'}, ...
+DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff');
+semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1);
+DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M);
+semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1);
+
+legend({'PSK with Viterbi-Viterbi', ...
+        'DEPSK with Viterbi-Viterbi', ...
+        'DPSK', ...
+        'Theoretical PSK over AWGN', ...
+        'Theoretical DEPSK over AWGN', ...
+        'Theoretical DPSK over AWGN'}, ...
        'Location', 'southwest');
 
-
 title({'QPSK with phase nosie and correction', ...
-       strcat(num2str(numSymbs * log2(M) / 1e3), '~kbit, LO~', ...
-              num2str(linewidthLO / 1e6), '~MHz, Av~', num2str(avgSa), ...
-              '~Sa')});
+       strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
+              '}$~bits, LO~', ...
+              num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
+              num2str(avgSa), '~Sa')});
 grid on;
 xlabel('$E_b/N_0$ (dB)');
 ylabel('BER');