--- /dev/null
+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.
--- /dev/null
+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;
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;
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;
ylabel('BER');
formatFigure;
- %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
- % '.svg'));
-
- %scatterplot(rxFilt);
- %eyediagram(rxFilt, sps);
-
end
--- /dev/null
+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.
--- /dev/null
+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');
--- /dev/null
+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');
-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);
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
[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'));
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
-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);
+++ /dev/null
-%% 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
--- /dev/null
+function y = normalizeEnergy(x, numSymbs, e)
+ energy = sum(abs(x) .^ 2) / numSymbs;
+ y = x * sqrt(e / energy);
+end
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)';
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;
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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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;
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / 4yp.git / commitdiff