phasenoise1signal.m

   1 numSymbs = 10000;
   2 M = 4;
   3
   4 Rsym = 2.5e10; % symbol rate (sym/sec)
   5 rolloff = 0.25;
   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 modData = pskmod(data, M, 0, 'gray');
  28
  29 x = txFilter(modData, rolloff, span, sps);
  30
  31 linewidthTx = 0;%1e5; % Hz
  32 linewidthLO = 1e6; % Hz
  33 %%linewidthTx = Rsym * 1e-4; % Hz
  34 %%linewidthLO = Rsym * 1e-3; % Hz
  35
  36 [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp);
  37
  38 snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
  39 noiseEnergy = 10 ^ (-snr / 10);
  40
  41 y = awgn(xPN, snr, 'measured');
  42
  43 r = rxFilter(y, rolloff, span, sps);
  44
  45 [rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps);
  46 rPhaseEq = normalizeEnergy(rPhaseEq, numSymbs, 1 + noiseEnergy);
  47
  48 rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps);
  49 demodData = pskdemod(rSampled, M, 0, 'gray')';
  50
  51 [bitErrors, ber] = biterr(data, demodData)
  52
  53
  54
  55 figure(2);
  56 plot(-phiests);
  57 hold on;
  58 plot(pTxLO);
  59 legend('estimate', 'actual');
  60 title('Phase noise estimation');
  61 hold off;
  62
  63 figure(3);
  64 plot(t(1:length(x)), real(x));
  65 hold on;
  66 plot(t, real(rPhaseEq), 'r');
  67 %%sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  68 %%plot(sampledTimes, real(rSampled), 'x');
  69 legend('original signal', 'corrected received signal');
  70 title('Phase noise correction, linewidth 1 MHz, E_b/N_0=8 dB');
  71 ylabel('Real part of signals');
  72 axis([t(1) t(300) -Inf +Inf]);
  73 hold off;
  74 formatFigure;