| 1 | numSymbs = 5e4; |
| 2 | M = 4; |
| 3 | rolloff = 0.5; |
| 4 | |
| 5 | Rsym = 2.5e10; % symbol rate (sym/sec) |
| 6 | |
| 7 | span = 6; % filter span |
| 8 | sps = 8; % samples per symbol |
| 9 | |
| 10 | fs = Rsym * sps; % sampling freq (Hz) |
| 11 | Tsamp = 1 / fs; |
| 12 | |
| 13 | t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; |
| 14 | |
| 15 | |
| 16 | EbN0_db = 0:0.2:14; |
| 17 | EbN0 = 10 .^ (EbN0_db ./ 10); |
| 18 | |
| 19 | Es = 1; |
| 20 | Eb = Es / log2(M); |
| 21 | N0 = Eb ./ EbN0; |
| 22 | |
| 23 | EsN0 = EbN0 .* log2(M); |
| 24 | EsN0_db = 10 .* log10(EsN0); |
| 25 | |
| 26 | plotlen = length(EbN0); |
| 27 | |
| 28 | berPSK = zeros(1, plotlen); |
| 29 | berDEPSK = zeros(1, plotlen); |
| 30 | berDPSK = zeros(1, plotlen); |
| 31 | |
| 32 | data = randi([0 M - 1], numSymbs, 1); |
| 33 | |
| 34 | pskSym = pskmod(data, M, pi / M, 'gray'); |
| 35 | %% DEPSK: Part VII, M.G. Taylor (2009) |
| 36 | depskSym = pskmod(data, M, 0, 'gray'); |
| 37 | for i = 2:numSymbs |
| 38 | depskSym(i) = depskSym(i) * depskSym(i-1); |
| 39 | end |
| 40 | |
| 41 | dpskSym = dpskmod(data, M, pi / M, 'gray'); |
| 42 | |
| 43 | xPSK = txFilter(pskSym, rolloff, span, sps); |
| 44 | xDEPSK = txFilter(depskSym, rolloff, span, sps); |
| 45 | xDPSK = txFilter(dpskSym, rolloff, span, sps); |
| 46 | |
| 47 | linewidthTx = 0; % Hz |
| 48 | linewidthLO = 5e6; % Hz |
| 49 | %linewidthLO = Rsym * 1e-3; |
| 50 | |
| 51 | iterations = 1; |
| 52 | avgSa = 40; |
| 53 | |
| 54 | TsampOrig = Tsamp; |
| 55 | |
| 56 | for it = 1 : iterations |
| 57 | [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp); |
| 58 | [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp); |
| 59 | [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp); |
| 60 | |
| 61 | for i = 1:plotlen |
| 62 | Tsamp = TsampOrig; |
| 63 | sps = 8; |
| 64 | |
| 65 | snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); |
| 66 | noiseEnergy = 10 ^ (-snr / 10); |
| 67 | |
| 68 | yPSK = awgn(xPSKpn, snr, 'measured'); |
| 69 | yDEPSK = awgn(xDEPSKpn, snr, 'measured'); |
| 70 | yDPSK = awgn(xDPSKpn, snr, 'measured'); |
| 71 | |
| 72 | rPSK = rxFilter(yPSK, rolloff, span, sps); |
| 73 | rDEPSK = rxFilter(yDEPSK, rolloff, span, sps); |
| 74 | rDPSK = rxFilter(yDPSK, rolloff, span, sps); |
| 75 | |
| 76 | sps = 2; |
| 77 | Tsamp = TsampOrig * 4; |
| 78 | |
| 79 | rPSKSamp = rPSK(1:2:end); |
| 80 | rDEPSKSamp = rDEPSK(1:2:end); |
| 81 | rDPSKSamp = rDPSK(1:2:end); |
| 82 | |
| 83 | [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa); |
| 84 | [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa); |
| 85 | |
| 86 | demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').'; |
| 87 | %% The decoding method described in Taylor (2009) |
| 88 | %% works on the complex symbols, i.e. after taking |
| 89 | %% the nearest symbol in the constellation, but before |
| 90 | %% converting them back to integers/bits. |
| 91 | %% MATLAB's pskdemod() does not provide this intermediate |
| 92 | %% result, so to be lazy, a pskmod() call is performed |
| 93 | %% to obtain the complex symbols. |
| 94 | demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').'; |
| 95 | remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray'); |
| 96 | delayed = [1; remodDEPSK(1:end-1)]; |
| 97 | demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray'); |
| 98 | |
| 99 | demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray'); |
| 100 | |
| 101 | [~, ber] = biterr(data, demodPSK); |
| 102 | berPSK(i) = berPSK(i) + ber / iterations; |
| 103 | [~, ber] = biterr(data, demodDEPSK); |
| 104 | berDEPSK(i) = berDEPSK(i) + ber / iterations; |
| 105 | [~, ber] = biterr(data, demodDPSK); |
| 106 | berDPSK(i) = berDPSK(i) + ber / iterations; |
| 107 | |
| 108 | if EbN0_db(i) == 8 && it == 1 |
| 109 | figure(1234); |
| 110 | plot(repelem(-phiestsPSK, 8)); |
| 111 | hold on; |
| 112 | plot(pTxLoPSK); |
| 113 | legend('estimate', 'actual'); |
| 114 | hold off; |
| 115 | |
| 116 | figure(1); |
| 117 | scatterplot(rPSKSampEq); |
| 118 | title('rPSKSampEq'); |
| 119 | end |
| 120 | |
| 121 | end |
| 122 | end |
| 123 | |
| 124 | |
| 125 | figure(1); |
| 126 | clf; |
| 127 | |
| 128 | %% Plot simulated results |
| 129 | semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5); |
| 130 | hold on; |
| 131 | semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2); |
| 132 | semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5); |
| 133 | |
| 134 | theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); |
| 135 | DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff'); |
| 136 | semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1); |
| 137 | DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M); |
| 138 | semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1); |
| 139 | |
| 140 | legend({'PSK with Viterbi-Viterbi', ... |
| 141 | 'DEPSK with Viterbi-Viterbi', ... |
| 142 | 'DPSK', ... |
| 143 | 'Theoretical PSK over AWGN', ... |
| 144 | 'Theoretical DEPSK over AWGN', ... |
| 145 | 'Theoretical DPSK over AWGN'}, ... |
| 146 | 'Location', 'southwest'); |
| 147 | |
| 148 | title({'QPSK with phase nosie and correction', ... |
| 149 | strcat('$10^{', num2str(log10(numSymbs * log2(M))), ... |
| 150 | '}$~bits, LO~', ... |
| 151 | num2str(linewidthLO / 1e6), '~MHz, blocksize~', ... |
| 152 | num2str(avgSa), '~Sa')}); |
| 153 | grid on; |
| 154 | xlabel('$E_b/N_0$ (dB)'); |
| 155 | ylabel('BER'); |
| 156 | |
| 157 | formatFigure; |