5 Rsym = 2.5e10; % symbol rate (sym/sec)
7 span = 6; % filter span
8 sps = 8; % samples per symbol
10 fs = Rsym * sps; % sampling freq (Hz)
13 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
17 EbN0 = 10 .^ (EbN0_db ./ 10);
23 EsN0 = EbN0 .* log2(M);
24 EsN0_db = 10 .* log10(EsN0);
26 plotlen = length(EbN0);
28 berPSK = zeros(1, plotlen);
29 berDEPSK = zeros(1, plotlen);
30 berDPSK = zeros(1, plotlen);
32 data = randi([0 M - 1], numSymbs, 1);
34 pskSym = pskmod(data, M, pi / M, 'gray');
35 %% DEPSK: Part VII, M.G. Taylor (2009)
36 depskSym = pskmod(data, M, 0, 'gray');
38 depskSym(i) = depskSym(i) * depskSym(i-1);
41 dpskSym = dpskmod(data, M, pi / M, 'gray');
43 xPSK = txFilter(pskSym, rolloff, span, sps);
44 xDEPSK = txFilter(depskSym, rolloff, span, sps);
45 xDPSK = txFilter(dpskSym, rolloff, span, sps);
48 linewidthLO = 5e6; % Hz
49 %linewidthLO = Rsym * 1e-3;
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);
65 snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
66 noiseEnergy = 10 ^ (-snr / 10);
68 yPSK = awgn(xPSKpn, snr, 'measured');
69 yDEPSK = awgn(xDEPSKpn, snr, 'measured');
70 yDPSK = awgn(xDPSKpn, snr, 'measured');
72 rPSK = rxFilter(yPSK, rolloff, span, sps);
73 rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
74 rDPSK = rxFilter(yDPSK, rolloff, span, sps);
77 Tsamp = TsampOrig * 4;
79 rPSKSamp = rPSK(1:2:end);
80 rDEPSKSamp = rDEPSK(1:2:end);
81 rDPSKSamp = rDPSK(1:2:end);
83 [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
84 [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
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');
99 demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');
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;
108 if EbN0_db(i) == 8 && it == 1
110 plot(repelem(-phiestsPSK, 8));
113 legend('estimate', 'actual');
117 scatterplot(rPSKSampEq);
128 %% Plot simulated results
129 semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
131 semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
132 semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);
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);
140 legend({'PSK with Viterbi-Viterbi', ...
141 'DEPSK with Viterbi-Viterbi', ...
143 'Theoretical PSK over AWGN', ...
144 'Theoretical DEPSK over AWGN', ...
145 'Theoretical DPSK over AWGN'}, ...
146 'Location', 'southwest');
148 title({'QPSK with phase nosie and correction', ...
149 strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
151 num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
152 num2str(avgSa), '~Sa')});
154 xlabel('$E_b/N_0$ (dB)');