phasenoise1signal.m

   1 numSymbs = 10000;
   2 M = 4;
   3
   4 Rsym = 2.5e10; % symbol rate (sym/sec)
   5 rolloff = 0.5;
   6 span = 6; % filter span
   7 sps = 4; % samples per symbol
   8
   9 fs = Rsym * sps; % sampling freq (Hz)
  10 Tsamp = 1 / fs;
  11
  12 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
  13
  14
  15 EbN0_db = 8;
  16 EbN0 = 10 .^ (EbN0_db ./ 10);
  17
  18 Es = 1;
  19 Eb = Es / log2(M);
  20 N0 = Eb ./ EbN0;
  21
  22 EsN0 = EbN0 .* log2(M);
  23 EsN0_db = 10 .* log10(EsN0);
  24
  25
  26 data = randi([0 M - 1], numSymbs, 1);
  27 pskSym = pskmod(data, M, pi/M, 'gray');
  28 dpskSym = dpskmod(data, M, pi/M, 'gray');
  29
  30 xPSK = txFilter(pskSym, rolloff, span, sps);
  31 xDPSK = txFilter(dpskSym, rolloff, span, sps);
  32
  33 linewidthTx = 0; % Hz
  34 linewidthLO = 10e6; % Hz
  35
  36 [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
  37 [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
  38
  39
  40 snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
  41
  42 %%y = awgn(xPN, snr, 'measured');
  43 yPSK = awgn(xPSKpn, snr, 'measured');
  44 yDPSK = xDPSKpn;
  45
  46 rPSK = rxFilter(yPSK, rolloff, span, sps);
  47 rDPSK = rxFilter(yDPSK, rolloff, span, sps);
  48
  49
  50 rPSKSa = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  51 rDPSKSa = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  52
  53 [rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 40);
  54
  55 demodPSK = pskdemod(rPSKSaPhEq, M, pi/M, 'gray');
  56 demodDPSK = dpskdemod(rDPSKSa, M, pi/M, 'gray');
  57
  58
  59 [bitErrors, ber] = biterr(data, demodPSK.')
  60
  61
  62 figure(2);
  63 plot(t(1:40000), repelem(-phiestsPSK, sps));
  64 hold on;
  65 plot(t(1:40000), pTxLoPSK(1:40000));
  66 legend('estimate', 'actual');
  67 title('Phase noise estimation');
  68 hold off;
  69 return
  70 figure(3);
  71 plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)));
  72 hold on;
  73 %%plot(t, real(rPhaseEq), 'r');
  74 sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  75 %%plot(sampledTimes, real(rSampled), 'x');
  76 plot(sampledTimes, real(normalizeEnergy(rSaPhEq, numSymbs, 1)), 'x');
  77 legend('original signal', 'corrected received signal');
  78 title('Phase noise correction, linewidth 1 MHz, $E_b/N_0=8$ dB');
  79 ylabel('Real part of signals');
  80 axis([t(1) t(300) -Inf +Inf]);
  81 hold off;
  82 formatFigure;