%蒙特卡羅仿真（AWGN信道，硬判決，2PSK）
%
Eb=2;     %每比特的能量,能量越大，卷積譯碼的誤碼率越小
k=1;
G=[1 0 1;1 1 1];
input=[1 0 1 0 1 1 1 0];
dsource=cnv_encd(G,k,input);
N=length(dsource);
n=size(G,1);    %得到碼字長度
Ec=k*Eb/n;      %得到碼元能量

% FSK調制
Tb=1;
f1=1000/Tb;
f2=f1+1/Tb;
%f2=1500/Tb;
phi=pi/4;
NN=500;
t=0:Tb/(NN-1):Tb;
u1=cos(2*pi*f1*t);   %若發送0，則調制信號為u1
u2=cos(2*pi*f2*t);   %若發送1，則調制信號為u2

% 通過Rayleigh信道
chan = rayleighchan(1/1000,0) 

for j=1:N
    if(dsource(j)==0)
        r=cos(2*pi*f1*t+phi);
        y= filter(chan,r);     %衰落
    else
        r=cos(2*pi*f2*t+phi);
        y= filter(chan,r);     %衰落
    end;
    %計算相關輸出
           v1=sin(2*pi*f1*t);
           v2=sin(2*pi*f2*t);
           r1c(1)=y(1)*u1(1);
           r1s(1)=y(1)*v1(1);
           r2c(1)=y(1)*u2(1);
           r2s(1)=y(1)*v2(1);
             for k=2:NN
              r1c(k)=r1c(k-1)+y(k)*u1(k);
              r1s(k)=r1s(k-1)+y(k)*v1(k);
              r2c(k)=r2c(k-1)+y(k)*u2(k);
              r2s(k)=r2s(k-1)+y(k)*v2(k);
             end;
     %判決變量 
     r1=r1c(NN)^2+r1s(NN)^2;
     r2=r2c(NN)^2+r2s(NN)^2; 
     if(abs(r1)>abs(r2))      %判決 
         decis(j)=0;
     else
         decis(j)=1;
     end;
end;
%硬判決維特比譯碼
[decoder_output,survivor_state,cumulated_metric]=viterbi(G,k,decis);