msk.m

介绍了MSK调制解调仿真程序

%整个程序
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*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=1:n
      if (v(i)>0),sa1(i)=-1;
      else sa1(i)=1;
      end;
  end; 
  
    figure
  plot(sa)
  figure
  plot(sa1)
  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