5 Rsym = 2.5e10; % symbol rate (sym/sec)
7 span = 6; % filter span
8 sps = 4; % 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;
54 for it = 1 : iterations
55 [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
56 [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
57 [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
60 snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
61 noiseEnergy = 10 ^ (-snr / 10);
63 yPSK = awgn(xPSKpn, snr, 'measured');
64 yDEPSK = awgn(xDEPSKpn, snr, 'measured');
65 yDPSK = awgn(xDPSKpn, snr, 'measured');
67 rPSK = rxFilter(yPSK, rolloff, span, sps);
68 rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
69 rDPSK = rxFilter(yDPSK, rolloff, span, sps);
71 rPSKSamp = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
72 rDEPSKSamp = rDEPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
73 rDPSKSamp = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
75 [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
76 [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
78 demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').';
79 %% The decoding method described in Taylor (2009)
80 %% works on the complex symbols, i.e. after taking
81 %% the nearest symbol in the constellation, but before
82 %% converting them back to integers/bits.
83 %% MATLAB's pskdemod() does not provide this intermediate
84 %% result, so to be lazy, a pskmod() call is performed
85 %% to obtain the complex symbols.
86 demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').';
87 remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray');
88 delayed = [1; remodDEPSK(1:end-1)];
89 demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray');
91 demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');
93 [~, ber] = biterr(data, demodPSK);
94 berPSK(i) = berPSK(i) + ber / iterations;
95 [~, ber] = biterr(data, demodDEPSK);
96 berDEPSK(i) = berDEPSK(i) + ber / iterations;
97 [~, ber] = biterr(data, demodDPSK);
98 berDPSK(i) = berDPSK(i) + ber / iterations;
100 if EbN0_db(i) == 8 && it == 1
102 plot(repelem(-phiestsPSK, sps));
105 legend('estimate', 'actual');
109 scatterplot(rPSKSampEq);
120 %% Plot simulated results
121 semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
123 semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
124 semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);
126 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
127 DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff');
128 semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1);
129 DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M);
130 semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1);
132 legend({'PSK with Viterbi-Viterbi', ...
133 'DEPSK with Viterbi-Viterbi', ...
135 'Theoretical PSK over AWGN', ...
136 'Theoretical DEPSK over AWGN', ...
137 'Theoretical DPSK over AWGN'}, ...
138 'Location', 'southwest');
140 title({'QPSK with phase nosie and correction', ...
141 strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
143 num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
144 num2str(avgSa), '~Sa')});
146 xlabel('$E_b/N_0$ (dB)');