clear all;
close all;

N=100;
N_Sample=8;
Ts=1;
dt=Ts/N_Sample;
t=0:dt:(N*N_Sample-1)*dt;
gt=ones(1,N_Sample);
d=sign(randn(1,N));
a=NumExpand(d,N_Sample);
st=conv(a,gt);

ht1=gt;
rt1=conv(st,ht1);

ht2=5*sinc(5*(t-5)/Ts);
rt2=conv(st,ht2);

figure(1)
subplot(321)
stem(d);
axis([0 20 -2 2]);
xlabel('輸入數字序列');
subplot(322)
plot(t,st(1:length(t)));
axis([0 20 -2 2]);
xlabel('輸入雙極性NRZ數字序列');
subplot(323)
plot(t,[0 rt1(1:length(t)-1)/8])
axis([0 20 -2 2]);
xlabel('方波濾波后輸出');
subplot(324);
dd=rt1(N_Sample:N_Sample:end);
ddd=NumExpand(dd,N_Sample);
stem(t,ddd(1:length(t))/8);
axis([0 20 -2 2]);
xlabel('方波濾波后抽樣輸出');
subplot(325)
plot(t-5,[0 rt2(1:length(t)-1)]/8);
axis([0 20 -2 2]);
xlabel('t/Ts   理想低通濾波后輸出');
subplot(326)
dd=rt2(N_Sample-1:N_Sample:end);
ddd=NumExpand(dd,N_Sample);
stem(t-5,ddd(1:length(t))/8);
axis([0 20 -2 2]);
xlabel('t/Ts   理想低通濾波后抽樣輸出');

%誤碼率性能仿真
EbN0dB=0:0.5:10;
N0=10.^(-EbN0dB/10);
sigma=sqrt(N0/2);
Pb=0.5*erfc(sqrt(1./N0));

for n=1:length(EbN0dB)
    a=sign(randn(1,100000));
    rk=a+sigma(n)*randn(1,100000);
    dec_a=sign(rk);
    ber(n)=sum(abs(a-dec_a)/2)/length(a);
end
figure(2)
semilogy(EbN0dB,Pb);
hold;
semilogy(EbN0dB,ber,'rd-');
legend('理論值','仿真結果');
title('誤碼率性能仿真');
xlabel('Eb/N0(dB)');ylabel('Pb');

%眼圖仿真
Ts=1;
N_Sample=17;
Eye_Num=7;
alpha=1;
N_data=1000;
dt=Ts/N_Sample;
t=-3*Ts:dt:3*Ts;

aa=sign(randn(1,N_data));
aaa=NumExpand(aa,N_Sample);
htt=sinc(t/Ts).*(cos(alpha*pi*t/Ts))./(1-4*alpha^2*t.^2/Ts^2+eps);
st=conv(aaa,htt);
st=awgn(st,30);
tt=-3*Ts:dt:(N_data+3)*N_Sample*dt-dt;
figure(3)
subplot(211);
plot(tt,st);
axis([0 20 -1.2 1.2]);
xlabel('t/Ts');title('基帶信號');
subplot(212)
ss=zeros(1,Eye_Num*N_Sample);
ttt=0:dt:(Eye_Num*N_Sample-1)*dt;
for k=3:50
    ss=st(k*N_Sample+1:(k+Eye_Num)*N_Sample);
    drawnow;
    plot(ttt,ss);
    hold on;
end
xlabel('t/Ts SNR=30dB');title('基帶信號眼圖');