function idpamsys(N_length,noisestd)
%idpamsys.m:  idealized baseband pam
%March 2008; Adapted from idsys.m from Johnson & Sethares
% N_length is number of 4-PAM symbols in simulation
% noisestd = wideband AWGN standard deviation

%clear all
clf

%TRANSMITTER
% Generate random 4-PAM (+/-1, +/-3) sequence
%N_length=1000;
m=m_pam(N_length,4); % Generate 4-PAM of length N_length samples
N=length(m);
% zero pad T-spaced symbol sequence to create upsampled T/M-spaced
% sequence of scaled T-spaced pulses (with T = 1 time unit)
M=100; mup=zeros(1,N*M); mup(1:M:end)=m; % oversampling factor
% Hamming pulse filter with T/M-spaced impulse response
p=hamming(M)/sqrt(pow(hamming(M))); % blip pulse of width M, unit norm
x=filter(p,1,mup);                  % convolve pulse shape with data

figure(1), plotspec(x,1/M)          % baseband signal spectrum
title('Transmit Signal Spectrum')
subplot(2,1,1)
title('Transmit PAM Signal')
% AWGN CHANNEL
%noisestd=1.0;
x1=x+sqrt(M)*noisestd*randn(1,N*M); % Add wide-band AWGN samples
%RECEIVER
figure(2), plotspec(x1,1/M) % spectrum of noisy received signal
title('Spectrum of Noisy Received Signal With Wideband AWGN')
subplot(2,1,1)
title('Noisy Received Signal')
% extract upsampled pulses using correlation implemented as a convolving filter
y=filter(fliplr(p)/M,1,x1); % matched filter; normalized
figure(5) % Compare transmit signal, noisy received signal, matched filter output
subplot(3,1,1)
plot(x)
title('Transmitted PAM Signal')
subplot(3,1,2)
plot(x1)
title('Noisy Received Signal')
subplot(3,1,3)
plot(y)
title('Output of Matched Filter')
figure(3), plotspec(y,1/M)
title('Spectrum of Matched Filter Output')
% set delay to one symbol = M samples, and increment by M
z=y(M:M:end); %+noisestd*randn(1,N);   % downsample to symbol rate; add noise
figure(4), plot([1:length(z)],z,'.',[1:length(m)],m,'x') % soft decisions
title('Transmitted symbols, x, Noisy Received Samples, \cdot')
% decision device and symbol matching performance assessment
mprime=quantalph(z,[-3,-1,1,3])';    % quantize to +/-1 and +/-3 alphabet
cluster_variance=(mprime-z)*(mprime-z)'/length(mprime), % cluster variance
lmp=length(mprime);
percentage_symbol_errors=100*sum(abs(sign(mprime-m(1:lmp))))/lmp % symb err
snr_dB=10*log10(5/noisestd^2), p_error_ml=1.5*0.5*erfc(sqrt(5/(2*5*noisestd^2)))


