n=10000;
count=100;
N=n/count;
r=rand(1,N-1);
for i=1:N-1
    if (r(i)>=0.5)
      r(i)=1;
   else r(i)=-1;
      end;
end;
for i=1:n-count
    R(i)=r(((i-1)-mod((i-1),count))/count+1);
end;
a(1)=1;
for i=2:N
    if r(i-1)==1
     a(i)=-a(i-1);
     else a(i)=a(i-1);
      end;
end;
for i=1:n
  A(i)=a(((i-1)-mod((i-1),count))/count+1);
end;
for i=1:N
    if(mod(i,2))
    l(i)=a(i);
    l(i+1)=a(i);
    else q(i)=a(i);
     q(i-1)=a(i);
     end;
 end;
for i=1:n
 I(i)=l(((i-1)-mod((i-1),count))/count+1);
 end;
for i=1:n 
  Q(i)=q(((i-1)-mod((i-1),count))/count+1);
end;
ts=0.0001;
Ts=ts*count;
fs=1/ts;
fc=5/(2*Ts);
t0=Ts*N-ts;
f=1/(4*Ts);
df=0.3;
t=[0:ts:t0];
tI=[-Ts:ts:t0-Ts];
tQ=[0:ts:t0];
tQ_R=[0:ts:t0-Ts];
c=2*pi*f;
u=I.*cos(c*tI);
v=Q.*sin(c*tQ);
U=u.*cos(2*pi*fc*tI);
V=v.*sin(2*pi*fc*tQ);
for i=count+1:n
  W(i)=U(i)+V(i-count);
end;
figure(1)
subplot(4,1,1);
plot(tQ_R,R,'k');
axis([-Ts,0.1,-2,2]);
title('数字信号');
grid on
subplot(4,1,2);
plot(tI,A,'k');
axis([-Ts,0.1,-2,2]);
title('差分编码后的信号');
grid on
subplot(4,1,3);
plot(tI,I,'k');
axis([-Ts,0.1,-2,2]);
title('同相调制信号I');
grid on
subplot(4,1,4);
 plot(tQ,Q,'k');
 axis([-Ts,0.1,-2,2]);
 title('正交调制信号Q');
 grid on
 figure(2)
 subplot(3,1,1);
 plot(tI,U,'k');
 axis([-Ts,0.1,-2,2]);
 title('Icos(wct)cos(wt)');
 grid on
 subplot(3,1,2);
 plot(tI,V,'k');
 axis([-Ts,0.1,-2,2]);
 title('Qsin(wct)sin(wt)');
 grid on
 subplot(3,1,3);
  plot(tI,W,'k');
  axis([-Ts,0.1,-2,2]);
  title('MSK信号');
  grid on
  noise=0.15*randn(1,n);
  s=ones(1,length(cos(2*pi*fc*t)));
  S=W+noise;
  [S,s,df]=fftseq(s,ts,df);
  S=S/fs;
  f=[0:df1:df1*(length(s)-1)]-fs/2;
  n_cutoff=100;
  H=zeros(size(f));
  for i=floor((fc-n_cutoff)/df1):floor((fc+n_cutoff)/df1)
      H(i)=1;
  end;
  for i=(65536-floor((fc-n_cutoff)/df1)):-1:(65536-floor((fc+n_cutoff)/df1)),H(i)=1;
  end;
  DEM=H.*S;
  dem=real(ifft(DEM))*fs;
  dm=dem(1:length(tQ));
  ds1=dm.*cos(2*pi*fc*tQ+pi*tQ/(2*Ts));
  ds2=dm.*cos(2*pi*fc*tQ-pi*tQ/(2*Ts));
  [DW1,ww,df1]=fftseq(ds1,ts,df);
  [DW2,ww,df1]=fftseq(ds2,ts,df);
  DW1=DW1/fs;
  DW2=DW2/fs;
  n_cutoff1=floor(200/df1);
  f1=[0:df1:df1*(length(ww)-1)]-fs/2;
  H1=zeros(size(f1));
  H1(1:n_cutoff1)=2*ones(1,n_cutoff1);
  H1(length(f1)-n_cutoff1+1:length(f1))=2*ones(1,n_cutoff1);
  DEM1=H1.*DW1;
  DEM2=H1.*DW2;
  s1=real(ifft(DEM1))*fs;
  s11=s1(1:length(tQ));
  s2=real(ifft(DEM2))*fs;
  s22=s2(1:length(tQ));
  for i=1:n
      if (s11(i)>=0.95)||(s11(i)<=-0.95),sa(i)=-1;
          else sa(i)=1;
       end;
end;
for i=1:n
    if (s22(i)>=0.95)||(s22(i)<=-0.95),sb(i)=1;
      else sb(i)=-1;
    end;
end;
 for i=0:N-1
     u1(i+1)=sa(i*count+count/2);
     u2(i+1)=sb(i*count+count/2);
end;
for j=1:N
    if (mod(j,2)==0) out(j)=u2(j);
    else out(j)=u1(j);
    end;
end;
for i=1:n
    OUT(i)=out(((i-1)-mod((i-1),count))/count+1);
end;
figure(3)
subplot(3,1,1);
plot(tQ,s,'k');
axis([0,0.1,-2,2]);
title(' 加噪声后的信号S');
axis([0,0.1,-2,2]);
grid on
subplot(3,1,2);
 plot(tQ,dm,'k');
 title(' 进过带通后的信号dm');
 axis([0,0.1,-2,2]);
 grid on
 subplot(3,1,3);
  plot(tQ,ds1,'k');
  title('?乘以载波后的信号ds1');
  axis([0,0.1,-2,2]);
 grid on 
 figure(4)
 subplot(3,1,1);
 plot(tQ,ds2,'k');
 title(' 乘以载波后的信号ds2');
 axis([0,0.1,-2,2]);
 grid on
 subplot(3,1,2);
 plot(tQ,s11,'k');
 title(' 经过低通后的信号s11');
 axis([0,0.1,-2,2]);
 grid on
 subplot(3,1,3);
 plot(tQ,s22,'k');
 title(' 经过低通后的信号s22');
 axis([0,0.1,-2,2]);
 grid on
 figure(5)
 subplot(3,1,1);
 plot(f,fftshift(abs(H)),'k');
 title('带通滤波器H');
 grid on
 subplot(3,1,2);
 plot(f1,fftshift(abs(H1)),'k');
 title('低通滤波器H1');
 grid on
 subplot(3,1,3);
 plot(tQ,sa,'k');
  title(' 经过低通后的信号sa');
  axis([0,0.1,-2,2]); 
  grid on
  figur(6)
  subplot(3,1,1);
  plot(tQ,sb,'k');
  title(' 经过低通后的信号sb');
  axis([0,0.1,-2,2]);
  grid on
  subplot(3,1,2);
  plot(tQ_R,R,'k');
  axis([0,1,-2,2]);
  title('调制信号');
  grid on
  subplot(3,1,3);
  plot(tQ_R,OUT(count+1:n),'k');
  title('解调信号');
  axis([0,1,-2,2]);
  grid on
  count1=0;
for i=1:N-1
   if(r(i)~=out(i+1))
     count1=count1+1;
    end;
end;
Pe=count1/(N-1)*100;
sprintf('MSK错误概率为：％g％％',Pe)
     d