From: Adrian Iain Lam Date: Fri, 2 Nov 2018 01:15:42 +0000 (+0000) Subject: Chromatic dispersion and line width phase noise X-Git-Url: https://adrianiainlam.tk/git/?p=4yp.git;a=commitdiff_plain;h=1eeb62fbc496ed5c170d199143ad53e28122d29c 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 --- 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