%用到的自編函數：快速傅立葉變換 
function[M,m,df]=fftseq(m,ts,df) 
fs=1/ts; 
n1=fs/df; 
n2=length(m); 
n=2^(max(nextpow2(n1),nextpow2(n2))); 
M=fft(m,n); 
m=[m,zeros(1,n-n2)]; 
df=fs/n; 
 
%整個程序 
n=10000; 
count=100; 
N=n/count;   %N=100 
r=rand(1,N-1); % 產生一個行向量，共9999個 
 
for i=1:N-1, 
    if (r(i)>=0.5) 
        r(i)=1; 
    else r(i)=-1; 
    end; 
end;           %將這9999個數整數化為 1，-1 
 
for i=1:n-count          % 取9900個 
    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; 
 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 
 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 
  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(tQ,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=W+noise; 
  [S1,s,df1]=fftseq(S,ts,df);   
  S1=S1/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=S1.*H; 
  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 
  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 
  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 
  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 
  figure 
  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) 
 
 
 
%頻率間隔的圖 
  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; 
r([1:10])=[1 -1 -1 1 1 1 -1 1 -1 -1]; 
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=3/(2*Ts); 
 t0=Ts*N-ts; 
 f=1/(4*Ts); 
 df=0.3; 
 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 
 subplot(3,1,1); 
 plot(tQ_R,R,'k');axis([0,0.1,-2,2]); 
 %title('數字信號'); 
  grid on  
  subplot(3,1,2); 
   plot(tI,W,'k');axis([0,0.1,-2,2]); 
  %title('MSK信號'); 
grid on