From 8449f934a70109b11c63373ca459dce3ac854cec Mon Sep 17 00:00:00 2001 From: Adrian Iain Lam Date: Fri, 26 Oct 2018 06:17:11 +0100 Subject: [PATCH] Added baseband and passband simulations. Passband is buggy. --- RRC_PSK_BER_SNR.m | 94 --------------------------------- baseband.m | 115 +++++++++++++++++++++++++++++++++++++++++ passband.m | 151 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 266 insertions(+), 94 deletions(-) delete mode 100644 RRC_PSK_BER_SNR.m create mode 100644 baseband.m create mode 100644 passband.m diff --git a/RRC_PSK_BER_SNR.m b/RRC_PSK_BER_SNR.m deleted file mode 100644 index 73d67f3..0000000 --- a/RRC_PSK_BER_SNR.m +++ /dev/null @@ -1,94 +0,0 @@ -function RRC_PSK_BER_SNR(rolloff, M, numSymbs) - %% Set defaults for inputs - if nargin < 3 - numSymbs = 1000; - end - if nargin < 2 - M = 2; - end - if nargin < 1 - rolloff = 0.5; - end - - if isOctave() - pkg load communications - end - - %% https://www.mathworks.com/help/signal/ref/rcosdesign.html - %% https://www.mathworks.com/help/comm/ug/pulse-shaping-using-a-raised-cosine-filter.html - span = 6; % filter span - sps = 4; - - rrcFilter = rcosdesign(rolloff, span, sps, 'sqrt'); - - EbN0_db = 0:0.2:10; - 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'); - - txSig = upfirdn(modData, rrcFilter, sps); - - for i = 1:plotlen - snr = EbN0_db(i) + 10 * log10(log2(M));% - 10 * log10(sps); % why sps? - rxSig = awgn(txSig, snr); - - rxFilt = upfirdn(rxSig, rrcFilter, 1, sps); - rxFilt = rxFilt(span + 1 : end - span); % remove filter delay - - demodData = pskdemod(rxFilt, M, 0, 'gray'); - - [bitErrors, ber(i)] = biterr(data, demodData); - end - - fig1 = 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 - - title(strcat(num2str(M), '-PSK RRC with Gray code')); - grid on; - xlabel('$E_b/N_0$ (dB)'); - ylabel('BER'); - - formatFigure; - %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ... - % '.svg')); - - %scatterplot(rxFilt); - %eyediagram(rxFilt, sps); - -end diff --git a/baseband.m b/baseband.m new file mode 100644 index 0000000..faf37f9 --- /dev/null +++ b/baseband.m @@ -0,0 +1,115 @@ +function baseband(rolloff, M, numSymbs) + %% Set defaults for inputs + if nargin < 3 + numSymbs = 1000; + end + if nargin < 2 + M = 2; + end + if nargin < 1 + 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; + 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'); + + xBaseband = txFilter([modData; zeros(span, 1)]); + + + + for i = 1:plotlen + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why 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); + + 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); + 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 + + title(strcat(num2str(M), '-PSK with Gray code')); + grid on; + xlabel('$E_b/N_0$ (dB)'); + ylabel('BER'); + + formatFigure; + %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ... + % '.svg')); + + %scatterplot(rxFilt); + %eyediagram(rxFilt, sps); + +end diff --git a/passband.m b/passband.m new file mode 100644 index 0000000..7abd749 --- /dev/null +++ b/passband.m @@ -0,0 +1,151 @@ +function passband(rolloff, M, numSymbs) + %% Set defaults for inputs + if nargin < 3 + numSymbs = 1000; + end + if nargin < 2 + M = 2; + end + if nargin < 1 + 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; + 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'); + + xBaseband = txFilter([modData; zeros(span, 1)]); + + %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]; + + carrier = zeros(length(t), 1); + for i = 1 : length(txLOFreq) + carrier = carrier + ... + sqrt(2 * txLOEnergy(i)) * exp(j * 2 * pi * txLOFreq(i) * t); + end + + xPassband = normalizeEnergy... + (real(xBaseband .* carrier(1:length(xBaseband))), numSymbs, 1); + + sum(abs(xPassband) .^ 2) / numSymbs + input('pause') + + + for i = 1:plotlen + snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why sps? + noiseEnergy = 10 ^ (-snr / 10); + + + yPassband = awgn(xPassband, snr, 'measured'); + + + rBaseband = rxFilter([yPassband .* 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); + + + 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); + 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 + + title(strcat(num2str(M), '-PSK RRC with Gray code')); + grid on; + xlabel('$E_b/N_0$ (dB)'); + 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 -- 2.7.4