仿真程序.txt

信号调制识别：在SNR=10db

function main()
global fc fs  fm fb  N  t  snr  E  T 
N=4000;  snr=10;  E=10; T=11;
fc=10000;fs=80000;
t=0:1/fs:(N-1)/fs; 
fm=500;ft=cos(2*pi*fm*t);
fb=2000;
typeset=['am  ';'dsb ';'lsb ';'usb ';'fm  ';'2ask';'4ask';'2fsk';'4fsk';'2psk';'4psk'];
confusion_matrix=zeros(T,T);
for signal=1:T
   for times=1:E
       type=deblank(typeset(signal,:));
       y=signal_produced(type,N);
       yy=signal_received(y,snr);
       char=parameter(yy);
       decision=decide(char);
       confusion_matrix(signal,decision)=confusion_matrix(signal,decision)+1;
   end
end
fprintf('SNR=%ddB,Confusion Matrix:\n',snr);
for r=1:T
    for c=1:T
        fprintf('%6d',confusion_matrix(r,c));
    end
    fprintf('\n');
end
*********************************************************************************************
function y=signal_produced(x,N)   
% Set parameters. 
fc=10000;fs=80000;
t=0:1/fs:(N-1)/fs; 
fm=500;ft=cos(2*pi*fm*t);
fb=2000;               
% Modulated signals.
switch  x
    
   case 'cw'
         y=cos(2*pi*fc*t);             %  signal of CW
   case  'am'
          ka=0.2;                    % index of AM
          y=(1+ka*ft).*cos(2*pi*fc*t);  %  signal Of AM 
   case 'dsb'
         y=ft.*cos(2*pi*fc*t);         % signal of DSB
   case 'lsb'
         y= ft.*cos(2*pi*fc*t)+imag(hilbert(ft)).*sin(2*pi*fc*t);   %  signal of SSB
   case 'usb'
         y= ft.*cos(2*pi*fc*t)-imag(hilbert(ft)).*sin(2*pi*fc*t);   %  signal of USB 
   case 'fm'
         kf=50;                      % index of FM
         y=cos(2*pi*fc*t+kf*sin(2*pi*fm*t));     %  signal of FM 
   case '2ask'
         n1=fs/fb;                    % The number of sample per symbol
         L=N/n1;               % The number of symbols
         PN=randsrc(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);
             end
         end
         for j=1:N
             if a(j)==1
                b(j)=1;
             else
                b(j)=0;
             end
         end
         y=b.*cos(2*pi*fc*t);    %  signal of 2ASK 
     case '4ask'
         n1=fs/fb;                    % The number of sample per symbol
         L=N/n1;               % The number of symbols
         PN=rand(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);
             end
         end
         for j=1:N
             if a(j)<0.25
                b(j)=0;
             elseif  a(j)<0.5
                b(j)=0.25;
             elseif  a(j)<0.75
                b(j)=0.50;
             else  b(j)=0.75;
             end
         end
         y=b.*cos(2*pi*fc*t);    %  signal of 4ASK
      case '2fsk'      
         n1=fs/fb;                    % The number of sample per symbol
         L=N/n1;               % The number of symbols
         PN=randsrc(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);
             end
         end
         f12=5000;
         y=cos(2*pi*fc*t+f12*2*pi*a.*t);            %  signal of 2FSK
     case '4fsk'      
         n1=fs/fb;                    % The number of sample per symbol
         L=N/n1;               % The number of symbols
         PN=randsrc(1,L);PN1=randsrc(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);u(n1*k+i)=PN1(k+1);
             end
         end
         f12=5000;f34=2500;                             
         y=cos(2*pi*fc*t+f12*2*pi*a.*t+f34*2*pi*u.*t);            %  signal of 4FSK..
      case '2psk'
        n1=fs/fb;                    % The number of sample per symbol
        L=N/n1;               % The number of symbols
        PN=randsrc(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);
             end
         end
         y=a.*cos(2*pi*fc*t);    %  signal of 2PSK
     case '4psk'
         n1=fs/fb;                    % The number of sample per symbol
         L=N/n1;               % The number of symbols
         PN=rand(1,L);
         for k=0:L-1
             for i=1:n1
                 a(n1*k+i)=PN(k+1);
             end
         end
         for j=1:N
             if a(j)<0.25
                b(j)=0;
             elseif  a(j)<0.5
                b(j)=1;
             elseif  a(j)<0.75
                b(j)=2;
             else  b(j)=3;
             end
         end
         y=cos(2*pi*fc*t+b.*pi/2);     %  signal of 4PSk.
         % case '4psk'
         %n1=fs/fb;                    % The number of sample per symbol
         %L=N/n1;               % The number of symbols
          %PN=randsrc(1,L);PN1=randsrc(1,L);
         %for k=0:L-1
             %for i=1:n1
                 %a(n1*k+i)=PN(k+1);b(n1*k+i)=PN1(k+1);
                 %end
                 %end
         %y=a.*cos(2*pi*fc*t)+b.*sin(2*pi*fc*t);     %  signal of 4PSk.
         
 end
*****************************************************************************************
function signal_figure1()
N=400;T=11;
fc=10000;fs=80000;
t=0:1/fs:(N-1)/fs; 
fm=500;ft=cos(2*pi*fm*t);
fb=2000;
typeset=['am  ';'dsb ';'lsb ';'usb ';'fm  ';'2ask';'4ask';'2fsk';'4fsk';'2psk';'4psk'];
for signal=1:T
       type=deblank(typeset(signal,:));
       y=signal_produced(type,N);
       %df1=fs/(length(y)-1);            
       %f=[0:df1:df1*(length(y)-1)]-fs/2; %???
       %f=[0:N-1]*fs;  %???
       %f=0:(N-1);  %???
       sig_spec=abs(fft(y));
       z=hilbert(y);   z1=hilbert(y).*exp(-j*2*pi*fc*t); 
       a=abs(z);
       phi_nl=angle(z1); 
       %PHI_NL=fft(phi);
       %fre=ifft(-j*2*pi*f.*PHI_NL); %???
       fre=diff(phi_nl)/(2*pi/fs);
       for i=1:(N-1)
           t1(i)=t(i);
       end
       figure(signal);
       subplot(2,3,1);
       plot(t,y);
       xlabel('t/s');
       title([upper(type),' signal']);
       subplot(2,3,2);
       %plot(f,sig_spec);  %???
       plot(sig_spec);
       xlabel('f/Hz');
       title([upper(type),' spectrum']);
       subplot(2,3,3);
       plot(t,a);
       xlabel('t/s');
       title([upper(type),' instaneous amplititude']);
       subplot(2,3,4);
       plot(t,phi_nl);
       xlabel('t/s');
       title([upper(type),' instaneous phase']);
       subplot(2,3,5);
       plot(t1,fre);
       xlabel('t/s');
       title([upper(type),' instaneous frequency']);
 end
*****************************************************************************************
function char=parameter(yy)
N=4000; 
fc=10000;fs=80000;
t=0:1/fs:(N-1)/fs;
fb=2000;    % guess  carrier frequency  ???

%%%%%%%insaneous amplititude%%%%%%instaneous phase%%%%%%%%%instaneous frequency%%%% 
z=hilbert(yy);   z1=hilbert(yy).*exp(-j*2*pi*fc*t); 
a=abs(z);
phi=angle(z);                  

ma=mean(a);  an=a/ma; acn=an-1;  

%phi_nl=phi-mean(phi);  ???
%phi_uw=unwrap(phi);
%for i=1:N
%phi_nl(i)=phi_uw(i)-2*pi*fc*i/fs;   ???
%end
phi_nl=angle(z1);

%PHI_NL=fft(phi_nl);
%df1=fs/(length(y)-1);            
%f=[0:df1:df1*(length(y)-1)]-fs/2;    ???
%f=[0:N-1]*fs;       ???
%f=0:(N-1);          ???
%fre=ifft(-j*2*pi*f.*PHI_NL);    ???  
%for i=1:(N-1)
%fre(i)=fs/2/pi*(phi_nl(i+1)-phi_nl(i));  ???
%end
fre=diff(phi_nl)/(2*pi/fs);


mf=mean(fre);  fmi=fre-mf;  fn=fmi/fb;   

%%%%%%%%%%%%  non-weak  segment %%%%%%%%%%%%%%%%
at=1;  k_non=0;                              
for k=1:(N-1)
    if an(k)>at
        k_non=k_non+1;
        acn_non(k_non)=acn(k);
        phi_nl_non(k_non)=phi_nl(k);
        fn_non(k_non)=fn(k);
        %else phi=pi/2;   ???
     end                           
end
        
%%%%pl%%%%%%pu%%%%%%%
fcn=fc*N/fs-1;
 s=fft(yy);pl=0; pu=0;
 for i=1:fcn
     pl=pl+(abs(s(i)))^2;
     pu=pu+(abs(s(i+fcn+1)))^2;
 end
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%parameter%%%%%%%%%%%%%%%%%%%%
 char1=max((fft(acn)).^2)/N;
 char2=std(abs(phi_nl_non),1);
 char3=std(phi_nl_non,1);
 char4=(pl-pu)/(pl+pu);
 char5=std(abs(acn),1);
 char6=std(abs(fn_non),1);
 char7=std(acn_non,1);
 char8=mean(acn.^4)/(mean(acn.^2))^2;
 char9=mean(fn.^4)/(mean(fn.^2))^2;
 char=[char1 char2 char3 char4 char5 char6 char7 char8 char9];
**************************************************************************************
function decision=decide(char)
if char(3)<=pi/6
    if char(8)>2.15 
        decision=1;    % AM
    elseif char(5)>0.25
        decision=7;  % 4ASk
      else decision=6; % 2ASK
      end
  elseif abs(char(4))>0.6
      if char(4)>0 decision=3; % LSB
      else decision=4; % USB 
      end
     elseif char(2)<=pi/5.5           
          if char(7)>0.25 decision=2; % DSB      
          else decision=10; % 2PSK
          end
       elseif  char(1)>2.5 decision=11; % 4PSK
        elseif char(9)>2.03 decision=5; % FM 
          elseif  char(6)>0.4 decision=9; % 4FSK
            else decision=8; % 2FSK 
  end
*******************************************************************************************
function yy=signal_received(y,snr)
snr_lin=10^(snr/10);
signal_power=(norm(y)^2)/length(y);   % signal_power=var(y)  ???
noise_power=signal_power/snr_lin;
noise_std=sqrt(noise_power);
noise=noise_std*randn(1,length(y));
yy=y+noise;               % signal+noise
%  band limited ???