%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author : Krishna Pillai
% Email : krishna@dsplog.com
% Version : 1.0
% Date : 14 June 2009
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script for generting transmit waveforms in a minimum shift
% keying, passing through Additive White Gaussian Noise and
% employing coherent demodulation at the receiver and then
% computing the number of errors
clear
N = 5*10^5; % number of bits or symbols
fsHz = 1; % sampling period
T = 4; % symbol duration
Eb_N0_dB = [0:10]; % multiple Eb/N0 values
ct = cos(pi*[-T:N*T-1]/(2*T));
st = sin(pi*[-T:N*T-1]/(2*T));
for ii = 1:length(Eb_N0_dB)
% MSK Transmitter
ipBit = rand(1,N)>0.5; % generating 0,1 with equal probability
ipMod = 2*ipBit - 1; % BPSK modulation 0 -> -1, 1 -> 0
ai = kron(ipMod(1:2:end),ones(1,2*T)); % even bits
aq = kron(ipMod(2:2:end),ones(1,2*T)); % odd bits
ai = [ai zeros(1,T) ]; % padding with zero to make the matrix dimension match
aq = [zeros(1,T) aq ]; % adding delay of T for Q-arm
% MSK transmit waveform
xt = 1/sqrt(T)*[ai.*ct + j*aq.*st];
% Additive White Gaussian Noise
nt = 1/sqrt(2)*[randn(1,N*T+T) + j*randn(1,N*T+T)]; % white gaussian noise, 0dB variance
% Noise addition
yt = xt + 10^(-Eb_N0_dB(ii)/20)*nt; % additive white gaussian noise
%% MSK receiver
% multiplying with cosine and sine waveforms
xE = conv(real(yt).*ct,ones(1,2*T));
xO = conv(imag(yt).*st,ones(1,2*T));
bHat = zeros(1,N);
bHat(1:2:end) = xE(2*T+1:2*T:end-2*T) > 0 ; % even bits
bHat(2:2:end) = xO(3*T+1:2*T:end-T) > 0 ; % odd bits
% counting the errors
nErr(ii) = size(find([ipBit - bHat]),2);
end
simBer = nErr/N; % simulated ber
theoryBer = 0.5*erfc(sqrt(10.^(Eb_N0_dB/10))); % theoretical ber
% plot
close all
figure
semilogy(Eb_N0_dB,theoryBer,'bs-','LineWidth',2);
hold on
semilogy(Eb_N0_dB,simBer,'mx-','LineWidth',2);
axis([0 10 10^-5 0.5])
grid on
legend('theory - bpsk', 'simulation - msk');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for MSK modulation');