% Simulation of suboptimum bandpass filter/delay-and-multiply demodulator
% with integrate-and-dump detection for DPSK; input bandpass filter is
% Butterworth with selectable bandwidth-bit period product
%
Eb_N0_dB_max = input('Enter maximum Eb/N0 in dB ');
Eb_N0_dB_min = input('Enter minimum Eb/N0 in dB ');
samp_bit = input('Enter number of samples per bit used in simulation ');
n_order = input('Enter order of Butterworth detection filter ');
BWT_bit = input('Enter filter bandwidth normalized by bit rate ');
N_bits = input('Enter total number of bits in simulation ')
clf											%  Clear any plots left over from previous runs
ss = sign(rand(1,N_bits)-.5);				%	Generate sequence of random +-1s (a digit for each bit)
data = 0.5*(ss+1);							%	Logical data is sequence of 1s and 0s
data_diff_enc = diff_enc(data);			    %  Differentially encode data for DPSK
s = 2*data_diff_enc-1;						%  Generate bipolar data for modulation
T_bit = 1;									%	Arbitrarily take bit time as 1 second
BW = BWT_bit/T_bit;							%	Compute filter bandwidth from BW*T_bit
[num,den] = butter(n_order,2*BW/samp_bit);	%	Obtain filter num/den coefficients
Eb_N0_dB = [];									%  Ensure that plotting arrays are empty
Perror = [];
k = 0;
for kk = Eb_N0_dB_min:Eb_N0_dB_max		    %  Loop for simulating detector for each desired Eb/N0
    k = k+1;
    Eb_N0_dB(k) = kk;
    Eb_N0 = 10^(Eb_N0_dB(k)/10);			    %	Convert desired Eb/N0 from dB to ratio
	Eb = T_bit;								%	Bit energy is T_bit, assuming amplitude = 1
	N0 = Eb/Eb_N0;							%	Compute noise spectral density from Eb/N0
	del_t = T_bit/samp_bit;					%	Compute sampling interval
	sigma_n = sqrt(N0/(2*del_t));			%   Compute standard deviations of noise samples
	t=0:del_t:N_bits*T_bit-del_t;			%	Set up time axis
	sig = [];	sig = s(ones(samp_bit,1),:);			%  Build array whose columns are samp_bit long
	sig = sig(:);							%  Convert matrix where bit samples occupy columns to vector
	zi = [];								%  Make sure various arrays are empty
	bits_out = [];	y_det = [];	noise = sigma_n*randn(size(sig));	    %  Form sequence of Gaussian noise samples	[y,zf] = filter(num,den,sig+noise,zi);	%  Filter signal plus noise with chosen filter	zi = zf;								%  Save final values for future initial conditions in a block simulation
	y_ref = delay1(y,samp_bit);			    %  Reference signal for multiply is 1-bit delayed signal + noise
	y_mult = y.*y_ref;						%  Multiply received signal plus noise with 1 bit delayed version  
	bits_out = int_and_dump(y_mult,samp_bit,N_bits);	% Integrate-and-dump detector is a subprogram	error_array = abs(bits_out-data);	    %  Compare detected bits with actual input data stream; errors = 1s
	error_array(1:5) = 0;					%  Don't include first 5 bits due to transients	ss = sum(error_array);					%  Sum to get total number of errors
	Perror(k) = ss/(N_bits-5);				%  Subtract 5 from denominator to account for initial errors set = 0
end											%  End of Eb/N0 loop 
disp('E_b/N_0, dB;  P_E')					%  Display computed probability of error values versus Eb/N0
disp([Eb_N0_dB' Perror'])
%  Plot simulated bit error probabilities versus Eb/N0
semilogy(Eb_N0_dB, Perror,'--'), grid, xlabel('E_b/N_0; dB'), ylabel('P_E'), hold
title('Simulation of BEP for delay-and-multiply detector with Butterworth prefilter for DPSK')
%  Plot theoretical bit error probability for optimum DPSK detector
semilogy(Eb_N0_dB, 0.5*exp(-10.^(Eb_N0_dB/10)))
%  Plot approximate theoretical result for suboptimum detector
semilogy(Eb_N0_dB, qfn(sqrt(10.^(Eb_N0_dB/10))),'-.')legend(['Simulation; BT = ', num2str(BWT_bit),'; ', num2str(N_bits),' bits'],'Theory; optimum differential detector','Theory; delay/multiply detector', 3)