From 1eeb62fbc496ed5c170d199143ad53e28122d29c Mon Sep 17 00:00:00 2001 From: Adrian Iain Lam Date: Fri, 2 Nov 2018 01:15:42 +0000 Subject: [PATCH] Chromatic dispersion and line width phase noise Implemented CD and phase noise simulations and their corrections; files written to (1) inspect waveform behaviour with time, and (2) plot BER/SNR graphs. For CD, an extra file was written to demonstrate (check) the pulse-broadening behaviour of the power of the signal. Code cleanup: move some code out to separate files for reuse; renamed some variables for consistency; removed 3F4 PSK upper bounds in graphs (bound too loose, not useful); removed isOctave due to the current heavy dependence on functions not supported by Octave. TODO: put CD and phase noise together --- CDCompensation.m | 30 +++++++++++ CD_AWGN.m | 75 ++++++++++++++++++++++++++ baseband.m | 58 +++----------------- chromaticDispersion.m | 37 +++++++++++++ chromaticDispersion1Signal.m | 63 ++++++++++++++++++++++ chromaticDispersionTest.m | 23 ++++++++ discretePSK_BER_SNR.m | 38 ++++---------- formatFigure.m | 8 ++- isOctave.m | 15 ------ normalizeEnergy.m | 4 ++ passband.m | 88 ++++--------------------------- phaseNoise.m | 16 ++++++ phaseNoiseCorr.m | 23 ++++++++ phasenoise1signal.m | 74 ++++++++++++++++++++++++++ phasenoise_AWGN.m | 122 +++++++++++++++++++++++++++++++++++++++++++ rxFilter.m | 20 +++++++ theoreticalPSK.m | 20 +++++++ txFilter.m | 17 ++++++ 18 files changed, 553 insertions(+), 178 deletions(-) create mode 100644 CDCompensation.m create mode 100644 CD_AWGN.m create mode 100644 chromaticDispersion.m create mode 100644 chromaticDispersion1Signal.m create mode 100644 chromaticDispersionTest.m delete mode 100644 isOctave.m create mode 100644 normalizeEnergy.m create mode 100644 phaseNoise.m create mode 100644 phaseNoiseCorr.m create mode 100644 phasenoise1signal.m create mode 100644 phasenoise_AWGN.m create mode 100644 rxFilter.m create mode 100644 theoreticalPSK.m create mode 100644 txFilter.m diff --git a/CDCompensation.m b/CDCompensation.m new file mode 100644 index 0000000..0fec7c8 --- /dev/null +++ b/CDCompensation.m @@ -0,0 +1,30 @@ +function yCDComp = CDCompensation(y, D, lambda, z, Tsamp) + %% Chromatic dispersion compensation. + %% Params: + %% - y: received waveform with CD + %% - D: dispersion coefficient (ps / (nm km)) + %% - lambda: wavelength (nm) + %% - z: length of fibre (km) + %% - Tsamp: sampling time (s) + %% Output: + %% - yCDComp: y after performing CD compensation + + %% Convert everything to SI base units + c = 299792458; % m/s + D = D * 1e-6; % s/m^2 + lambda = lambda * 1e-9; % m + z = z * 1e3; % m + + %% Implementing Eq. (9) in [1]. + N = 2 * floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2)) + 1; + k = -floor(N / 2) : floor(N / 2); + + % h: FIR filter + h = sqrt(j * c * Tsamp^2 / (D * lambda^2 * z)) * ... + exp(-j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z)); + + yCDComp = upfirdn(y, h); + yCDComp = yCDComp(N:end); % truncate filter transients +end +%% References +%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008. diff --git a/CD_AWGN.m b/CD_AWGN.m new file mode 100644 index 0000000..1ea6ad8 --- /dev/null +++ b/CD_AWGN.m @@ -0,0 +1,75 @@ +numSymbs = 10000; +M = 4; + +Rsym = 2.5e10; % symbol rate (sym/sec) + +rolloff = 0.25; +span = 6; % filter span +sps = 4; % samples per symbol + +fs = Rsym * sps; % sampling freq (Hz) +Tsamp = 1 / fs; + +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; + +EbN0_db = 0:0.2:14; +EbN0 = 10 .^ (EbN0_db ./ 10); + +Es = 1; +Eb = Es / log2(M); +N0 = Eb ./ EbN0; + +EsN0 = EbN0 .* log2(M); +EsN0_db = 10 .* log10(EsN0); + +plotlen = length(EbN0); + +ber = zeros(1, plotlen); + +data = randi([0 M - 1], numSymbs, 1); +modData = pskmod(data, M, 0, 'gray'); +x = txFilter(modData, rolloff, span, sps); + +%% Simulate chromatic dispersion +D = 20; % ps / (nm km) +lambda = 1550; % nm +z = 10; % km + +xCD = chromaticDispersion(x, D, lambda, z, Tsamp); +xCD = normalizeEnergy(xCD, numSymbs, 1); + + +for i = 1:plotlen + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); + noiseEnergy = 10 ^ (-snr / 10); + + y = awgn(xCD, snr, 'measured'); + + yCDComp = CDCompensation(y, D, lambda, z, Tsamp); + + r = rxFilter(yCDComp, rolloff, span, sps); + %% normalize energy + %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy); + + rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps); + demodData = pskdemod(rSampled, M, 0, 'gray'); + + [bitErrors, ber(i)] = biterr(data, demodData); +end + +figure(1); +clf; + +%% Plot simulated results +semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); +hold on; + +theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); +legend({'CD + AWGN + CD compensation', 'AWGN only'}, 'Location', 'southwest'); + +title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation')); +grid on; +xlabel('$E_b/N_0$ (dB)'); +ylabel('BER'); + +formatFigure; diff --git a/baseband.m b/baseband.m index faf37f9..f3534fe 100644 --- a/baseband.m +++ b/baseband.m @@ -10,30 +10,16 @@ function baseband(rolloff, M, numSymbs) rolloff = 0.5; end - %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html Rsym = 1e6; % symbol rate (sym/sec) span = 6; % filter span sps = 4; % samples per symbol - txFilter = comm.RaisedCosineTransmitFilter... - ('Shape', 'Square root', ... - 'RolloffFactor', rolloff, ... - 'FilterSpanInSymbols', span, ... - 'OutputSamplesPerSymbol', sps); - rxFilter = comm.RaisedCosineReceiveFilter... - ('Shape', 'Square root', ... - 'RolloffFactor', rolloff, ... - 'FilterSpanInSymbols', span, ... - 'InputSamplesPerSymbol', sps, ... - 'DecimationFactor', 1); - fs = Rsym * sps; % sampling freq (Hz) t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; - EbN0_db = 0:0.2:10; EbN0 = 10 .^ (EbN0_db ./ 10); Es = 1; @@ -49,56 +35,30 @@ function baseband(rolloff, M, numSymbs) data = randi([0 M - 1], numSymbs, 1); modData = pskmod(data, M, 0, 'gray'); - xBaseband = txFilter([modData; zeros(span, 1)]); - - + xBaseband = txFilter(modData, rolloff, span, sps); for i = 1:plotlen - snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why sps? + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); noiseEnergy = 10 ^ (-snr / 10); yBaseband = awgn(xBaseband, snr, 'measured'); - rBaseband = rxFilter([yBaseband; zeros(span, 1)]); - %% truncate filter transients - rBaseband = rBaseband(span * sps / 2 + 1 : end); - %% normalize to unit energy - rBasebandEnergy = sum(abs(rBaseband) .^ 2) / numSymbs; - rBaseband = rBaseband .* sqrt((1 + noiseEnergy) / rBasebandEnergy); + rBaseband = rxFilter(yBaseband, rolloff, span, sps); rSampled = rBaseband(sps*span/2+1:sps:(numSymbs+span/2)*sps); - demodData = pskdemod(rSampled, M, 0, 'gray'); - [bitErrors, ber(i)] = biterr(data, demodData); end - fig1 = figure(1); + figure(1); clf; %% Plot simulated results semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); hold on; - %% Plot theoretical curve - %% BPSK: bit error when noise Nr > sqrt(Eb) - %% Pr(Nr > sqrt(Eb)) - %% = Pr(Z > sqrt(Eb) / sqrt(N0/2)) - %% - %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit - %% so BER is the same as BPSK (assuming Gray code) - if M == 2 || M == 4 - ber_th = qfunc(sqrt(2 * EbN0)); - semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1); - legend('Simulated', 'Discrete'); - else - %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary - %% Communication Systems using MATLAB (Equations - %% 7.3.18 and 7.3.19), Brooks/Cole. - ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M); - semilogy(EbN0_db, ber_ap, 'b', 'LineWidth', 1); - legend('Simulated', 'Discrete'); - end + theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); + legend('Simulated', 'Discrete'); title(strcat(num2str(M), '-PSK with Gray code')); grid on; @@ -106,10 +66,4 @@ function baseband(rolloff, M, numSymbs) ylabel('BER'); formatFigure; - %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ... - % '.svg')); - - %scatterplot(rxFilt); - %eyediagram(rxFilt, sps); - end diff --git a/chromaticDispersion.m b/chromaticDispersion.m new file mode 100644 index 0000000..49a49a3 --- /dev/null +++ b/chromaticDispersion.m @@ -0,0 +1,37 @@ +function [xCD, kstart] = chromaticDispersion(x, D, lambda, z, Tsamp) + %% Simulate chromatic dispersion. + %% Params: + %% - x: input waveform (pulse-shaped) + %% - D: dispersion coefficient (ps / (nm km)) + %% - lambda: wavelength (nm) + %% - z: length of fibre (km) + %% - Tsamp: sampling time (s) + %% Output: + %% - xCD: x after being dispersed. Energy of xCD is not normalized. + %% - kstart: starting index of the discrete signal + + %% Convert everything to SI base units + c = 299792458; % m/s + D = D * 1e-6; % s/m^2 + lambda = lambda * 1e-9; % m + z = z * 1e3; % m + + %% kmax: maximum k at which to sample the inpulse response. + %% See [1] Eq. (8), noting that the same omega (Eq. (7)) + %% can be used for dispersion and dispersion compensation, + %% as phi(t) is the same in both Eqs. (5) and (6). + kmax = floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2)); + k = -kmax : kmax; % index for discrete function + + t = k * Tsamp; + + % Impulse response + g = sqrt(c / (j * D * lambda^2 * z)) * ... + exp(j * pi * c / (D * lambda^2 * z) * t .^ 2); + + xCD = conv(x, g); + + kstart = 1 - kmax; +end +%% References +%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008. diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m new file mode 100644 index 0000000..db25965 --- /dev/null +++ b/chromaticDispersion1Signal.m @@ -0,0 +1,63 @@ +M = 4; +numSymbs = 1000; + +%% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html +Rsym = 2.5e10; % symbol rate (sym/sec) + +span = 6; % Tx/Rx filter span +rolloff = 0.25; % Tx/Rx RRC rolloff +sps = 4; % samples per symbol + + +fs = Rsym * sps; % sampling freq (Hz) +Tsamp = 1 / fs; + +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; + + +data = randi([0 M - 1], numSymbs, 1); +modData = pskmod(data, M, 0, 'gray'); +x = txFilter(modData, rolloff, span, sps); + +%% Simulate chromatic dispersion +D = 20; % ps / (nm km) +lambda = 1550; % nm +z = 1000; % km + +[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp); +xCD = normalizeEnergy(xCD, numSymbs, 1); + + +y = xCD; + + +yCDComp = CDCompensation(y, D, lambda, z, Tsamp); + +%% Compare original signal and compensated signal +figure(1); +subplot(211); +plot(real(x(1:300))); +hold on +plot(real(yCDComp(1:300))); +hold off +title('Real part'); +legend('original', 'dispersion compensated'); +subplot(212); +plot(imag(x(1:300))); +hold on +plot(imag(yCDComp(1:300))); +hold off +title('Imag part'); + + +r = rxFilter(yCDComp, rolloff, span, sps); +r = normalizeEnergy(r, numSymbs, 1); % Add noise energy if needed + +rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps); + +scatterplot(modData); +title('Constellation of original modulation'); +scatterplot(rSampled); +title('Constellation of sampled received waveform'); + +demodData = pskdemod(rSampled, M, 0, 'gray'); diff --git a/chromaticDispersionTest.m b/chromaticDispersionTest.m new file mode 100644 index 0000000..81af8b1 --- /dev/null +++ b/chromaticDispersionTest.m @@ -0,0 +1,23 @@ +Rsym = 2.5e10; % symbol/sec +sps = 4; % samples/symbol +fs = Rsym * sps; % sampling freq (Hz) +Tsamp = 1 / fs; + +x = 10 * normpdf(1:30, 10, 3); + +xCD = chromaticDispersion(x, 20, 1550, 30, Tsamp); + +figure(1); +subplot(211); +plot(abs(x) .^ 2); +axis([0 30 -Inf Inf]); +title('Before dispersion'); +ylabel('Power (arb. unit)'); +subplot(212); +plot(abs(xCD) .^ 2) +axis([20 50 -Inf Inf]); +title('After dispersion'); +xlabel('Time (arb. unit)'); +ylabel('Power (arb. unit)'); + +%saveas(gcf, 'dispersionTest.png'); diff --git a/discretePSK_BER_SNR.m b/discretePSK_BER_SNR.m index fd8849c..cdff01f 100644 --- a/discretePSK_BER_SNR.m +++ b/discretePSK_BER_SNR.m @@ -1,16 +1,12 @@ -function discretePSK_BER_SNR(M, numBits) +function discretePSK_BER_SNR(M, numSymbs) %% Set defaults for inputs if nargin < 2 - numBits = 1000; + numSymbs = 1000; end if nargin < 1 M = 2; end - if isOctave() - pkg load communications - end - EbN0_db = 0:0.2:10; EbN0 = 10 .^ (EbN0_db ./ 10); @@ -25,7 +21,7 @@ function discretePSK_BER_SNR(M, numBits) ber = zeros(1, plotlen); - data = randi([0, M - 1], 1, numBits); + data = randi([0, M - 1], 1, numSymbs); txsig = pskmod(data, M, 0, 'gray'); for i = 1:plotlen @@ -35,35 +31,19 @@ function discretePSK_BER_SNR(M, numBits) [bitErrors, ber(i)] = biterr(data, demodData); end - fig1 = figure(1); + figure(1); clf; %% Plot simulated results semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); hold on; - %% Plot theoretical curve - %% BPSK: bit error when noise Nr > sqrt(Eb) - %% Pr(Nr > sqrt(Eb)) - %% = Pr(Z > sqrt(Eb) / sqrt(N0/2)) - %% - %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit - %% so BER is the same as BPSK (assuming Gray code) + theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); + if M == 2 || M == 4 - ber_th = qfunc(sqrt(2 * EbN0)); - semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1); legend('Simulated', 'Theoretical'); else - %% Upper bound: R. Venkataramanan, Lent 2018, 3F4 Examples Paper 2 - %% (Question 5), CUED. - %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary - %% Communication Systems using MATLAB (Equations - %% 7.3.18 and 7.3.19), Brooks/Cole. - ber_ub = 2 * qfunc(sqrt(EbN0 * log2(M)) * sin(pi / M)); - ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M); - semilogy(EbN0_db, ber_ub, 'b', 'LineWidth', 1); - semilogy(EbN0_db, ber_ap, 'g', 'LineWidth', 1); - legend('Simulated', 'Upper bound', 'Approximation'); + legend('Simulated', 'Approximation'); end title(strcat(num2str(M), '-PSK with Gray code')); @@ -72,6 +52,6 @@ function discretePSK_BER_SNR(M, numBits) ylabel('BER'); formatFigure; - saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numBits), ... - '.svg')); + %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ... + % '.svg')); end diff --git a/formatFigure.m b/formatFigure.m index 7d5796a..a09586d 100644 --- a/formatFigure.m +++ b/formatFigure.m @@ -1,8 +1,6 @@ -if ~isOctave() - set(gcf, 'defaultTextInterpreter', 'Latex'); - set(gcf, 'defaultLegendInterpreter', 'Latex'); - set(gca, 'TickLabelInterpreter', 'Latex'); -end +set(gcf, 'defaultTextInterpreter', 'Latex'); +set(gcf, 'defaultLegendInterpreter', 'Latex'); +set(gca, 'TickLabelInterpreter', 'Latex'); set(gca, 'FontSize', 14); leg = findobj(gcf, 'Type', 'Legend'); set(leg, 'FontSize', 14); diff --git a/isOctave.m b/isOctave.m deleted file mode 100644 index 1c9c190..0000000 --- a/isOctave.m +++ /dev/null @@ -1,15 +0,0 @@ -%% Taken from the GNU Octave Manual, (c) 1996-2016 John W. Eaton -%% https://octave.org/doc/v4.0.3/How-to-distinguish-between-Octave-and-Matlab_003f.html - -%% -%% Return: true if the environment is Octave. -%% -function retval = isOctave - persistent cacheval; % speeds up repeated calls - - if isempty (cacheval) - cacheval = (exist ("OCTAVE_VERSION", "builtin") > 0); - end - - retval = cacheval; -end diff --git a/normalizeEnergy.m b/normalizeEnergy.m new file mode 100644 index 0000000..8c185dd --- /dev/null +++ b/normalizeEnergy.m @@ -0,0 +1,4 @@ +function y = normalizeEnergy(x, numSymbs, e) + energy = sum(abs(x) .^ 2) / numSymbs; + y = x * sqrt(e / energy); +end diff --git a/passband.m b/passband.m index 5424ea4..9e88af6 100644 --- a/passband.m +++ b/passband.m @@ -10,25 +10,12 @@ function passband(rolloff, M, numSymbs) rolloff = 0.5; end - %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html - Rsym = 1e6; % symbol rate (sym/sec) + Rsym = 2.5e10; % symbol rate (sym/sec) span = 6; % filter span sps = 4; % samples per symbol - txFilter = comm.RaisedCosineTransmitFilter... - ('Shape', 'Square root', ... - 'RolloffFactor', rolloff, ... - 'FilterSpanInSymbols', span, ... - 'OutputSamplesPerSymbol', sps); - rxFilter = comm.RaisedCosineReceiveFilter... - ('Shape', 'Square root', ... - 'RolloffFactor', rolloff, ... - 'FilterSpanInSymbols', span, ... - 'InputSamplesPerSymbol', sps, ... - 'DecimationFactor', 1); - fs = Rsym * sps; % sampling freq (Hz) t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; @@ -49,80 +36,39 @@ function passband(rolloff, M, numSymbs) ber = zeros(1, plotlen); - - data = randi([0 M - 1], numSymbs, 1); modData = pskmod(data, M, 0, 'gray'); - xBaseband = txFilter([modData; zeros(span, 1)]); - + xBaseband = txFilter(modData, rolloff, span, sps); - %fc = 2.5e6; % Carrier freq (Hz) - %carrier = sqrt(2) * exp(j * 2 * pi * fc * t); - - %xPassbandIdeal = normalizeEnergy... - % (real(xBaseband .* carrier(1:length(xBaseband))), numSymbs, 1); - - %txLOFreq = [2.49e6, 2.5e6, 2.51e6]; - %%txLOEnergy = [0.05, 0.9, 0.05]; - %txLOEnergy = [0 1 0]; - txLOFreq = [1e8]; - txLOEnergy = [1]; - - carrier = zeros(length(t), 1); - for i = 1 : length(txLOFreq) - carrier = carrier + ... - sqrt(2 * txLOEnergy(i)) * exp(j * 2 * pi * txLOFreq(i) * t); - end + fc = 1e12; % Carrier freq (Hz) + carrier = sqrt(2) * exp(j * 2 * pi * fc * t); - xPassband = normalizeEnergy... - (xBaseband .* carrier(1:length(xBaseband)), numSymbs, 1); + xPassband = xBaseband .* carrier(1:length(xBaseband)); for i = 1:plotlen - snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why sps? + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); noiseEnergy = 10 ^ (-snr / 10); yPassband = awgn(xPassband, snr, 'measured'); - - rBaseband = rxFilter([yPassband .* conj(carrier(1:length(yPassband))); zeros(span * sps, 1)]); - %% truncate filter transients - rBaseband = rBaseband(span * sps / 2 + 1 : end); - %% normalize energy - rBaseband = normalizeEnergy(rBaseband, numSymbs, 1 + noiseEnergy); + yLO = yPassband .* conj(carrier(1:length(yPassband))); + rBaseband = rxFilter(yLO, rolloff, span, sps); rSampled = rBaseband(sps*span/2+1:sps:(numSymbs + span/2) * sps); demodData = pskdemod(rSampled, M, 0, 'gray'); [bitErrors, ber(i)] = biterr(data, demodData); - end - fig1 = figure(1); + figure(1); clf; %% Plot simulated results semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); hold on; - %% Plot theoretical curve - %% BPSK: bit error when noise Nr > sqrt(Eb) - %% Pr(Nr > sqrt(Eb)) - %% = Pr(Z > sqrt(Eb) / sqrt(N0/2)) - %% - %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit - %% so BER is the same as BPSK (assuming Gray code) - if M == 2 || M == 4 - ber_th = qfunc(sqrt(2 * EbN0)); - semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1); - legend('Simulated RRC', 'Discrete'); - else - %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary - %% Communication Systems using MATLAB (Equations - %% 7.3.18 and 7.3.19), Brooks/Cole. - ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M); - semilogy(EbN0_db, ber_ap, 'b', 'LineWidth', 1); - legend('Simulated RRC', 'Discrete'); - end + theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); + legend('Simulated RRC', 'Discrete'); title(strcat(num2str(M), '-PSK RRC with Gray code')); grid on; @@ -130,16 +76,4 @@ function passband(rolloff, M, numSymbs) ylabel('BER'); formatFigure; - %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ... - % '.svg')); - - %scatterplot(rxFilt); - %eyediagram(rxFilt, sps); - -end - - -function y = normalizeEnergy(x, numSymbs, e) - energy = sum(abs(x) .^ 2) / numSymbs; - y = x * sqrt(e / energy); end diff --git a/phaseNoise.m b/phaseNoise.m new file mode 100644 index 0000000..82427b3 --- /dev/null +++ b/phaseNoise.m @@ -0,0 +1,16 @@ +function [xPN, phiTx_phiLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp) + %% fIF = 15e3; % fIF = fc - fLO + %% carrier = exp(-j * (2 * pi * fIF * t - phi); + %% Neglecting fIF for now... + + dphiTx = normrnd(0, sqrt(2 * pi * linewidthTx * Tsamp), length(x), 1); + dphiLO = normrnd(0, sqrt(2 * pi * linewidthLO * Tsamp), length(x), 1); + phiTx = cumsum(dphiTx); + phiLO = cumsum(dphiLO); + + phiTx_phiLO = phiTx - phiLO; + + pn = exp(-j * phiTx_phiLO); + + xPN = x .* pn; +end diff --git a/phaseNoiseCorr.m b/phaseNoiseCorr.m new file mode 100644 index 0000000..1f4ed45 --- /dev/null +++ b/phaseNoiseCorr.m @@ -0,0 +1,23 @@ +function [rPhaseEq, phiests] = phaseNoiseCorr(r, M, blocksize) + %% phase noise correction + phiests = zeros(1, length(r)); + rPhaseEq = zeros(1, length(r)); + for l = 1:blocksize:length(r) + block = r(l : min(l + blocksize - 1, length(r))); + + sum_M = sum(block .^ M); + phi_est = angle(sum_M) / M; % assume phase of 0 symbol is 0. + + if l > 1 + %% phase unwrapping + phi_prev = phiests(l - 1); + m = floor(0.5 + (phi_prev - phi_est) * M / (2 * pi)); + phi_est = phi_est + m * 2 * pi / M; + end + + phiests(l:min(l+blocksize-1, length(r))) = phi_est * ones(1, min(blocksize, length(r) - l+1)); + + block = block .* exp(j * -phi_est); + rPhaseEq(l : min(l + blocksize-1, length(r))) = block; + end +end diff --git a/phasenoise1signal.m b/phasenoise1signal.m new file mode 100644 index 0000000..e0ee39d --- /dev/null +++ b/phasenoise1signal.m @@ -0,0 +1,74 @@ +numSymbs = 10000; +M = 4; + +Rsym = 2.5e10; % symbol rate (sym/sec) +rolloff = 0.25; +span = 6; % filter span +sps = 4; % samples per symbol + +fs = Rsym * sps; % sampling freq (Hz) +Tsamp = 1 / fs; + +t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; + + +EbN0_db = 8; +EbN0 = 10 .^ (EbN0_db ./ 10); + +Es = 1; +Eb = Es / log2(M); +N0 = Eb ./ EbN0; + +EsN0 = EbN0 .* log2(M); +EsN0_db = 10 .* log10(EsN0); + + +data = randi([0 M - 1], numSymbs, 1); +modData = pskmod(data, M, 0, 'gray'); + +x = txFilter(modData, rolloff, span, sps); + +linewidthTx = 0;%1e5; % Hz +linewidthLO = 1e6; % Hz +%%linewidthTx = Rsym * 1e-4; % Hz +%%linewidthLO = Rsym * 1e-3; % Hz + +[xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp); + +snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps); +noiseEnergy = 10 ^ (-snr / 10); + +y = awgn(xPN, snr, 'measured'); + +r = rxFilter(y, rolloff, span, sps); + +[rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps); +rPhaseEq = normalizeEnergy(rPhaseEq, numSymbs, 1 + noiseEnergy); + +rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps); +demodData = pskdemod(rSampled, M, 0, 'gray')'; + +[bitErrors, ber] = biterr(data, demodData) + + + +figure(2); +plot(-phiests); +hold on; +plot(pTxLO); +legend('estimate', 'actual'); +title('Phase noise estimation'); +hold off; + +figure(3); +plot(t(1:length(x)), real(x)); +hold on; +plot(t, real(rPhaseEq), 'r'); +%%sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); +%%plot(sampledTimes, real(rSampled), 'x'); +legend('original signal', 'corrected received signal'); +title('Phase noise correction, linewidth 1 MHz, E_b/N_0=8 dB'); +ylabel('Real part of signals'); +axis([t(1) t(300) -Inf +Inf]); +hold off; +formatFigure; diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m new file mode 100644 index 0000000..a9c852a --- /dev/null +++ b/phasenoise_AWGN.m @@ -0,0 +1,122 @@ +function phasenoise_AWGN(rolloff, M, numSymbs) + %% Set defaults for inputs + if nargin < 3 + numSymbs = 1000; + end + if nargin < 2 + M = 2; + end + if nargin < 1 + rolloff = 0.25; + end + + plotted = 0; + + Rsym = 2.5e10; % symbol rate (sym/sec) + + span = 6; % filter span + sps = 4; % samples per symbol + + fs = Rsym * sps; % sampling freq (Hz) + Tsamp = 1 / fs; + + t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)'; + + + EbN0_db = 0:0.2:14; + EbN0 = 10 .^ (EbN0_db ./ 10); + + Es = 1; + Eb = Es / log2(M); + N0 = Eb ./ EbN0; + + EsN0 = EbN0 .* log2(M); + EsN0_db = 10 .* log10(EsN0); + + plotlen = length(EbN0); + + ber = zeros(1, plotlen); + + + data = randi([0 M - 1], numSymbs, 1); + modData = pskmod(data, M, 0, 'gray'); + + x = txFilter(modData, rolloff, span, sps); + + linewidthTx = 0;%1e5; % Hz + linewidthLO = 1e6; % Hz + %%linewidthTx = Rsym * 1e-4; % Hz + %%linewidthLO = Rsym * 1e-3; % Hz + + totalIterations = 1; + + for iter = 1:totalIterations + [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp); + + for i = 1:plotlen + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); + noiseEnergy = 10 ^ (-snr / 10); + + y = awgn(xPN, snr, 'measured'); + + r = rxFilter(y, rolloff, span, sps); + %% normalize energy + %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy); + + [rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps); + + rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps); + demodData = pskdemod(rSampled, M, 0, 'gray')'; + + %%[bitErrors, ber(i)] = biterr(data, demodData); + [zzz, thisBER] = biterr(data, demodData); + ber(i) = ber(i) + thisBER / totalIterations; + + + if plotted == 0 && EbN0_db(i) >6 && ber(i) > 1e-1 + plotted=1 + figure(1234); + plot(-phiests); + hold on; + plot(pTxLO); + legend('estimate', 'actual'); + hold off; + + figure(100); + %plot(t(1:length(x)), real(x)); + %%plot(t, real(x(1:length(t)))); + length(t) + length(x) + hold on + %%plot(t(1:length(xPhaseNoise)), real(xPhaseNoise)); + + plot(t, real(r), 'g') + + sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); + + plot(sampledTimes, real(rSampled), 'x') + hold off + + end + + end + end + + + figure(1); + clf; + + %% Plot simulated results + semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); + hold on; + + theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); + legend({'Simulated phase noise', 'Without phase noise'}, 'Location', 'southwest'); + + title(strcat(num2str(M), '-PSK with phase noise and correction')); + grid on; + xlabel('$E_b/N_0$ (dB)'); + ylabel('BER'); + + formatFigure; +end diff --git a/rxFilter.m b/rxFilter.m new file mode 100644 index 0000000..eb25a32 --- /dev/null +++ b/rxFilter.m @@ -0,0 +1,20 @@ +function r = rxFilter(y, rolloff, span, sps) + %% Receiver matched (root raised cosine) filter. + %% Inputs: + %% - y: received waveform + %% - rolloff: rolloff factor in root raised cosine filter. + %% - span: filter span (number of symbols) + %% - sps: samples per symbol + %% Output: + %% - r: filtered signal (energy not normalized), + %% normalize and then sample. + filter = comm.RaisedCosineReceiveFilter... + ('Shape', 'Square root', ... + 'RolloffFactor', rolloff, ... + 'FilterSpanInSymbols', span, ... + 'InputSamplesPerSymbol', sps, ... + 'DecimationFactor', 1); + + r = filter([y; zeros(span * sps, 1)]); + r = r(span * sps / 2 + 1 : end); % truncate filter transients +end diff --git a/theoreticalPSK.m b/theoreticalPSK.m new file mode 100644 index 0000000..95fbf97 --- /dev/null +++ b/theoreticalPSK.m @@ -0,0 +1,20 @@ +function theoreticalPSK(EbN0_db, M, varargin) + %% Plot theoretical curve + EbN0 = 10 .^ (EbN0_db ./ 10); + if M == 2 || M == 4 + %% BPSK: bit error when noise Nr > sqrt(Eb) + %% Pr(Nr > sqrt(Eb)) + %% = Pr(Z > sqrt(Eb) / sqrt(N0/2)) + %% + %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit + %% so BER is the same as BPSK (assuming Gray code) + ber_th = qfunc(sqrt(2 * EbN0)); + semilogy(EbN0_db, ber_th, varargin{:}); + else + %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary + %% Communication Systems using MATLAB (Equations + %% 7.3.18 and 7.3.19), Brooks/Cole. + ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M); + semilogy(EbN0_db, ber_ap, varargin{:}); + end +end diff --git a/txFilter.m b/txFilter.m new file mode 100644 index 0000000..49164d3 --- /dev/null +++ b/txFilter.m @@ -0,0 +1,17 @@ +function x = txFilter(modData, rolloff, span, sps) + %% Transmitter pulse-shaping (root raised cosine) filter. + %% Inputs: + %% - modData: modulated data + %% - rolloff: rolloff factor in root raised cosine filter. + %% - span: filter span (number of symbols) + %% - sps: samples per symbol + %% Output: + %% - x: pulse-shaped waveform + + filter = comm.RaisedCosineTransmitFilter... + ('Shape', 'Square root', ... + 'RolloffFactor', rolloff, ... + 'FilterSpanInSymbols', span, ... + 'OutputSamplesPerSymbol', sps); + x = filter([modData; zeros(span, 1)]); +end -- 2.7.4