function demoPAMrect(a,M,duty,dev)
% This function demonstrates the generation
% of a PAM signal from the symbol amplitudes
% in the vector a, using an upsampling factor
% of M. duty is the percent duty cycle, 0<duty<=100
% dev is the noise standard deviation.
% Sample run: demoPAMrect([3 3 -1 1 -3 1],100,50,.25)
aup=upsample(a,M); % Upsampled Data Symbol Sequence
%g=hamming(M,'symmetric')';
g=[ones(1,duty) zeros(1,M-duty)];
% g=hanning(M,'symmetric')'; % Transmit Pulse
%g=hanning(M,'periodic')'; % Transmit Pulse
scale=sqrt(sum(g.*g));
g=g/scale; % Normalize to unit energy
h=[0 fliplr(g)]; % Matched Filter h(n) = g(M-n)
s=conv(g,h); % Zero ISI Pulse
x=conv(g,aup); % Filter Data to Form Transmitted Signal
x1=filter(h,1,aup); % Alternative filtering command
y=conv(h,x); % Apply Matched Filter
y1=filter(fliplr(g)/(M*sum(g.^2)/length(g)),1,x1); % Alternative MF sum(x.^2)/length(x);
figure(1) % Plot g(n), h(n) and s(n)
subplot(2,1,1)
plot([0:length(g)-1],g,[0:length(h)-1],h,'--')
xlabel('Time Index, n')
ylabel('g(n), h(n)')
title('Transmit and Receive Filter Impulse Responses')
legend('Transmit Filter, g(n)','Receive Filter, h(n) = g(M-n)')
subplot(2,1,2)
plot(s)
xlabel('Time Index, n')
ylabel('s(n) = g(n)*h(n)')
title('Zero ISI Pulse')
figure(2) % Plot data symbols, transmit signal, and MF output
subplot(2,1,1)
stem([0:M*length(a)-1],[aup zeros(1,M*length(a)-length(aup))],'.')
xlabel('Time Index, n')
ylabel('a_{up}(n)')
title('Input Symbol Sequence')
subplot(2,1,2)
plot([0:length(x)-1],scale*x,'--',[0:length(y)-1],y) %,[0:length(y1)-1],y1)
xlabel('Time Index, n')
ylabel('x(n), y(n)')
title('Transmitted Signal and Matched Filter Output')
legend('Trans. Signal, x(n)','MF Output, y(n)')
hold on
stem([1:length(a)]*M,y(1+M:M:1+M*(length(a))),'.')
hold off
% Add noise
NumSim=20;
ahat=zeros(NumSim,length(a));
yout=zeros(NumSim,length(a));
figure(3)
for nn=1:20
z=dev*randn(1,length(x));
yn=conv(h,z+x);
ahat(nn,:)=uquantvec(yn(1+M:M:1+M*(length(a))),2,4);
yout(nn,:)=yn(1+M:M:1+M*(length(a)));
%subplot(2,1,1)
%stem([0:M*length(a)-1],[aup zeros(1,M*length(a)-length(aup))],'.')
%xlabel('Time Index, n')
%ylabel('a_{up}(n)')
%title('Input Symbol Sequence')
%subplot(2,1,2)
plot([0:length(y)-1],y,'--',[0:length(yn)-1],yn,0,5,0,-5) %,[0:length(y1)-1],y1)
xlabel('Time Index, n')
ylabel('y(n), y_{noisy}(n)')
title('Ideal and Noisy Matched Filter Outputs')
%legend('Ideal Output, y(n)','Noisy Output, y_{noisy}(n)')
legend('Ideal','Noisy')
hold on
stem([1:length(a)]*M,yn(1+M:M:1+M*(length(a))),'.')
hold off
pause(0.25)
end
figure(4)
% plot([1:length(a)]*M,a(1,:),'*')
% figure(5)
% plot([1:length(a)]*M,ahat(1,:),'*',[1:length(a)]*M,yout(1,:),'.')
hold on
for nn=1:NumSim
 %   plot([1:length(a)]*M,ahat(nn,:),'.')
    plot([1:length(a)]*M,yout(nn,:),'.')
end
hold off
xlabel('Time Index, n')
ylabel('Matched Filter Output Samples')
title('Ideal (*) and Noisy (\cdot) Matched Filter Output Samples')
end
function y = uquantvec(x,delta,num_levels)
% Uniform scalar quantizer with step size = delta
% num_levels is the number of quantization output
% levels, and the output levels are mid-points
% of the quantization regions
% if x is a vector, y is the quantized version of all samples in x.
if(num_levels < 1)
  error('num_levels must be positive integer')
end
N=length(x);
y=zeros(1,N);
for n=1:N
 if mod(num_levels,2) == 0
  s=1;
  if x(n)<0
    s=-1;
  end
  i=floor(abs(x(n))/delta);
  i=min(i,(num_levels/2)-1);
  y(n)=s*(i+0.5)*delta;
 else
  i=round(x(n)/delta);
  i=max(min(i,(num_levels-1)/2),-(num_levels-1)/2);
  y(n)=i*delta;
 end
end
end
    