Added adaptive CMA equalizer, "fixed" CD and phase noise

[4yp.git] / phasenoise_AWGN.m

diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m
index a9c852a..d5185be 100644 (file)
--- 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;