matlab.txt

这个题目是利用人工神经网络实现对数字调制信号的识别

clear all;
close all;
fc = 2000; %载波频率
fs = 40000; %采样速率
k=2;
code_size=15*round(k*fs/fc);             %  信息码元长度
t0=5.5;  %信号长度
%code_length=16O; %码元个数
Ns=256; %采样点个数
fd = 125;  %符号速率
ts=1/fs;  %采样周期
M=600 %码元个数
ti=1/fd; %码元间隔
N=ti/ts;
t=[0:ts:t0];
select=menu('调制方式','2ASK','2FSK','2PSK','4ASK','4FSK','4PSK');
switch select
case 1,   % 2ASK signal
x=randint(1,M);
m=sin(2*pi*fc*t);
y=ones(1,M*N);
for i=1:M
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end
end
T=zeros(6,50);
T(1,1:50)=1; 
case 2,%2FSK signal
x=randint(1,M);  
m1=sin(2*pi*fc*t);
m2=sin(2*pi*2*fc*t);
y=zeros(1,M*N);
for i=1:M
   if x(i)==1;
    for j=1:N;
        y((i-1)*N+j)=x(i)*m1(j);
    end
   elseif  x(i)==0;
     for j=1:N;
       y((i-1)*N+j)=(1-x(i))*m2(j);
     end 
  end
end
T=zeros(6,50);
T(2,1:50)=1; 
case 3,     %2PSK signal,
x=randint(1,M); 
m1=sin(2*pi*fc*t);
m2=sin(2*pi*fc*t+pi);
y=zeros(1,M*N);
for i=1:M
   if x(i)==1;
    for j=1:N;
        y((i-1)*N+j)=x(i)*m1(j);
   end
   elseif  x(i)==0;
     for j=1:N;
       y((i-1)*N+j)=(1-x(i))*m2(j);
     end 
   end
end
T=zeros(6,50);
T(3,1:50)=1;
case 4,      % 4ASK signal
x1=randint(1,M);
x2=randint(1,M);
x3=randint(1,M);
x=x1+x2+x3;
m=sin(2*pi*fc*t);
y=ones(1,M*N);
for i=1:M
    if x(i)==0;
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end
elseif x(i)==1;
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end
elseif x(i)==2;
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end
elseif x(i)==3;
    for j=1:N
       y((i-1)*N+j)=x(i)*m(j);
    end  
  end  
end
T=zeros(6,50);
T(4,1:50)=1; 
case 5,    % 4FSK signal
x1=randint(1,M);
x2=randint(1,M);
x3=randint(1,M);
x=x1+x2+x3;
m1=sin(2*pi*fc*t);
m2=sin(2*pi*2*fc*t);
m3=sin(2*pi*3*fc*t);
m4=sin(2*pi*4*fc*t);
y=zeros(1,M*N);
for i=1:M
    if x(i)==0;
    for j=1:N
       y((i-1)*N+j)=(1-x(i))*m1(j);
    end
elseif x(i)==1;
    for j=1:N
       y((i-1)*N+j)=x(i)*m2(j);
    end
elseif x(i)==2;
    for j=1:N
       y((i-1)*N+j)=(x(i)-1)*m3(j);
    end
elseif x(i)==3;
    for j=1:N
       y((i-1)*N+j)=(x(i)-2)*m4(j);
    end  
  end  
end
T=zeros(6,50);
T(5,1:50)=1;
case 6,    %4PSK signal
x1=randint(1,M);
x2=randint(1,M);
x3=randint(1,M);
x=x1+x2+x3;
m1=sin(2*pi*fc*t);
m2=sin(2*pi*fc*t+pi/2);
m3=sin(2*pi*fc*t+pi);
m4=sin(2*pi*fc*t+3*pi/2);
y=zeros(1,M*N);
for i=1:M
    if x(i)==0;
    for j=1:N
       y((i-1)*N+j)=(1-x(i))*m1(j);
    end
elseif x(i)==1;
    for j=1:N
       y((i-1)*N+j)=x(i)*m2(j);
    end
elseif x(i)==2;
    for j=1:N
       y((i-1)*N+j)=(x(i)-1)*m3(j);
    end
elseif x(i)==3;
    for j=1:N
       y((i-1)*N+j)=(x(i)-2)*m4(j);
    end  
  end  
end
T=zeros(6,50);
T(6,1:50)=1;
end
SNR =20;                                      %定义信噪比，单位DB 
sigpow=mean(abs(y).^2);  %power of input signal
noisefac=10^(-SNR/10);
noise=randn(1,size(y,2));
noise=noise*(sqrt(sigpow*noisefac)/sqrt(mean(noise.^2)));  %产生所需的高斯噪声
ynoise=noise+y;  %加噪后的信号
for n=1:1:50
    m=n*Ns;
    x=(n-1)*Ns;
for i=x+1:m;  %提取信号段
    y0(i)=ynoise(i);
end
Y=fft(y0);       %调制信号的傅立叶变换
y1=hilbert(y0);   %实信号的解析式 
z=abs(y0);       %实信号的瞬时幅度
phase=angle(y1); %实信号的瞬时相位
add=0;           %求Rmax
for i=x+1:m;
   add=add+z(i);
end
ma=add/Ns;     %瞬时幅度的平均值
y2=z./ma  ;    %幅度比,即为文献中的an(i)
y3=y2-1;       %归一化瞬时幅度
y4=max(abs(y3));
y2n=y3./y4;    % 即为文献中的acn(i)    
s=fft(y2n);   
R=abs(s);
Rmax=max((R)/Ns).^2; %零中心归一化瞬时幅度的谱密度的最大值
%Fcn=fc.*Ns/fs-1; %求取谱对称性P
%Pl=0;Pu=0;  
%for i=1:Fcn
    %Pl=Pl+abs(Y(i)).*abs(Y(i));
    %Pu=Pu+abs(Y(i+Fcn+1)).*abs(Y(i+Fcn+1));
    %end
%P=abs((Pl-Pu)/(Pl+Pu));  %谱对称性P
Xcn=0;
Ycn=0;
for i=x+1:m;
    Xcn=Xcn+y2n(i).*y2n(i);
    Ycn=Ycn+abs(y2n(i));
end
Xcnav=Xcn/Ns;
Ycnav=(Ycn/Ns).*(Ycn/Ns);
deltaaa=sqrt(Xcnav-Ycnav);  %零中心归一化瞬时幅度绝对值得标准偏差        
if  phase(2+x)-phase(1+x)>pi;  %修正相位序列
    Ck(1+x)=-2*pi;
elseif  phase(1+x)-phase(2+x)>pi;
    Ck(1+x)=2*pi;
else Ck(1+x)=0;
end
for i=x+2:m-1;
  if   phase(i+1)-phase(i)>pi; 
              Ck(i)=Ck(i-1)-2*pi;
      elseif  phase(i)-phase(i+1)>pi
              Ck(i)=Ck(i-1)+2*pi;
      else
              Ck(i)=Ck(i-1);
   end
end
if  -phase(m)>pi;
    Ck(m)=Ck(m-1)-2*pi;
elseif phase(m)>pi;
    Ck(m)=Ck(m-1)+2*pi;
else Ck(m)=Ck(m-1);
end
phase1=phase+Ck %去相位卷叠后的相位序列
phasenl=phase1-2*pi*fc*i/fs; %非线性相位
    at=1; %判决门限电平
    a=0;   %求取零中心非弱信号段瞬时相位非线性分量绝对值的标准偏差和零中心非弱信号段瞬时相位非线性分量的标准偏差
    b=0;
    d=0;
    c=0;
for i=x+1:m;
     if y2(i)>at
         c=c+1;
         phasesquare(i)=phasenl(i).*phasenl(i);
         a=a+phasesquare(i);
         phaseabs(i)=abs(phasenl(i));
         b=b+phaseabs(i);
         d=d+phasenl(i)        
     end
end
a1=a/c;
b1=(b/c).*(b/c);
d1=(d/c).*(d/c);
deltaap=sqrt(a1-b1); %零中心非弱信号段瞬时相位非线性分量绝对值的标准偏差
deltadp=sqrt(a1-d1);%零中心非弱信号段瞬时相位非线性分量的标准偏差
for i=x+1:m;
    if iat;
            c=c+1;
         freqNsquare(i)=freqN(i).*freqN(i);
         a=a+freqNsquare(i);
         b=b+freqN(i);
        end
end
a1=a/c;
b1=(b/c)^2;
deltaaf=sqrt(a1-b1); %零中心归一化非弱信号段瞬时频率绝对值得标准偏差
   
    Pi=rand(8,50);
    P0=rand(8,50);  
    %Pi(9*n-8)=P;
    Pi(6*n-5)=Rmax;
    Pi(6*n-4)=deltaap;
    Pi(6*n-3)=deltadp;
    Pi(6*n-2)=deltaaa;
    Pi(6*n-1)=deltaaf;
end
%采用BP网络
%NEWCF--生成一个新的前向神经网络
%TRAIN--对网络进行训练
% 定义训练样本
%Pi为输入矢量
%T为目标矢量
%创建一个新的前向神经网络　
net=newff(minmax(Pi),[15,6],{'tansig','purelin'},'traingda');
%当前输入层权值和阈值
inputWeights=net.IW{1,1};
inputWeights ;
inputbias=net.b{1};
inputbias 
%当前网络层权值和阈值
layerWeights=net.LW{2,1};
layerWeights; 
layerbias=net.b{2};
layerbias; 
%设置训练参数
net.trainParam.show=50;
net.trainParam.lr=1;
net.trainParam.mc=0.5;
net.trainParam.epochs=2000;
net.trainParam.goal=1e-3;
%调用TRAINGDM算法训练网络
[net,tr]=train(net,Pi,T);
%对网络进行仿真
A=sim(net,Pi);
A
E=T-A;

T1=zeros(6,50);
T1(1,1:50)=1;
E1=0;
for c=2:1:200
    Eout=(T1(c)-A(c))^2;  
    E1=Eout+E1 ;
end   
T2=zeros(6,50);
T2(2,1:50)=1;
E2=0;
for c=2:1:200
    Eout=(T2(c)-A(c))^2;
    E2=Eout+E2;
end   
T3=zeros(6,50);
T3(3,1:50)=1;
E3=0;
for c=2:1:200
    Eout=(T3(c)-A(c))^2;  
    E3=Eout+E3;
end   
T4=zeros(6,50);
T4(4,1:50)=1;
E4=0;
for c=2:1:200
    Eout=(T4(c)-A(c))^2;
    E4=Eout+E4;
end   
T5=zeros(6,50);
T5(5,1:50)=1;
E5=0;
for c=2:1:200
    Eout=(T5(c)-A(c))^2;
    E5=Eout+E5;
end   
T6=zeros(6,50);
T6(6,1:50)=1;
E6=0;
for c=2:1:200
    Eout=(T6(c)-A(c))^2;
    E6=Eout+E6;
end   
E0=0;
if (E1>E2)
    E0=E2;
else E0=E1;
end
if (E3
    E0=E3;
end
if (E4
    E0=E4;
end
if (E5
    E0=E5;
end
if (E6
    E0=E6;
end
E0;
if (E0==E1)
   type=menu('输入信号是','2ASK信号'); 
end
if (E0==E2)
   type=menu('输入信号是','2FSK信号'); 
end
if(E0==E3)
   type=menu('输入信号是','2PSK信号'); 
end
if(E0==E4)
   type=menu('输入信号是','4ASK信号'); 
end
if(E0==E5)
   type=menu('输入信号是','4FSK信号'); 
end
if(E0==E6)
   type=menu('输入信号是','4PSK信号'); 
end
%计算正确识别率sita
sita=0;
j=0;
for nn=1:1:50
   Ee=(E(6*nn-5)^2)+(E(6*nn-4)^2)+(E(6*nn-3)^2)+(E(6*nn-2)^2)+(E(6*nn-1)^2)+(E(6*nn)^2);
   if Ee<0.01
       j=j+1
   end
end
sita=j/50;
type=menu('正确识别率为:',sita*100,'％');