--- /dev/null
+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