function [p]=syssim(snr_in_dB)
E=1;
SNR=exp(snr_in_dB*log(10)/10); %信噪比
sgma=E/sqrt(2*SNR);	 %噪聲的標準偏差
M=10000;  %碼元數 
L=32;  %每個碼元采樣點數
N=M*L;  %采樣點數
Rb=2;   %碼速率
Ts=1/Rb;  %碼元間隔
dt=Ts/L;  %時域采樣間隔
df=1/(N*dt);  %頻域采樣間隔
d=16;  %半個碼元采樣點數
T=N*dt;  %截短時寬
Bs=N*df/2;  %頻帶寬度
t=[-T/2+dt/2:dt:T/2];   %時域坐標
t1=[-T/2+dt/2:32*dt:T/2];   %時域坐標
f=[-Bs+df/2:df:Bs];   %頻域坐標

alpha=0.5;	%滾降系數
g1=sin(pi*t/Ts)./(pi*t/Ts);
g2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
g=g1.*g2;
G=sqrt(t2f(g,dt));  %發送和接收濾波器頻譜
Cf=1-0.5*exp(-2*i*pi*f*d);  %	多徑信道

Gf=G.^2.*Cf;  %系統總的頻譜
%%  產生序列
a = sign(randn(1,M));
b = sign(randn(1,M));
%%  a路經過發送端的成型濾波器和信道
imp=zeros(1,N);
imp(L/2:L:N)=a/dt;
S=t2f(imp,dt).*Gf;
s=f2t(S,dt);

s1=zeros(1,N);       %此段程序將經傅立葉變換得到的數字序列
s1(N/2+1:N)=s(1:N/2);%在中軸反轉,如111-1-11,經變換變成-1-11111
s1(1:N/2)=s(N/2+1:N);%因此應把它反過來,這主要是由
s1=real(s1);         %f2t和t2f函數引起的.
gr=s1+sgma*randn(1,N);  %信號加噪聲，實現SNR
gr1=gr(16:32:N);  %取每個碼元中間的一個樣本
%%  b路經過發送端的成型濾波器和信道
imp=zeros(1,N);
imp(L/2:L:N)=b/dt;
S=t2f(imp,dt).*Gf;
s=f2t(S,dt);

s1=zeros(1,N);       %此段程序將經傅立葉變換得到的數字序列
s1(N/2+1:N)=s(1:N/2);%在中軸反轉,如111-1-11,經變換變成-1-11111
s1(1:N/2)=s(N/2+1:N);%因此應把它反過來,這主要是由
s1=real(s1);         %f2t和t2f函數引起的.
gr=s1+sgma*randn(1,N);  %信號加噪聲，實現SNR
gr2=gr(16:32:N);  %取每個碼元中間的一個樣本

%%  判決和誤碼概率的計算
numoferr=0;
for ii=1:M,
	%隨后為判決
    if (gr1(ii)<0),
		decis=-1;	 %判決為“0”
	else 
		decis=1;		%判決為“1”
	end;
	if(decis~=a(ii)),
		numoferr=numoferr+1;  %如果出錯，計數器加1
	end;
	if (gr2(ii)<0),
		decis=-1;	 %判決為“0”
	else 
		decis=1;		%判決為“1”
	end;
	if(decis~=b(ii)),
		numoferr=numoferr+1;  %如果出錯，計數器加1
	end;
end;
p=numoferr/2/M;		  %誤差概率估算