信噪比的定义很多，我们利用的是local signal-to-noise ratio，这个需要取信号谱线下一个小的频带内噪声功率，与直接除以噪声方差是不一样的。

   下面是一段计算饱和非线性系统的输出信噪比程序，不过解释不多，如果想看的很明白，可以email我，fabing.duan@gmail.com

Matlab code:

clear all

f=0.01;%signal frequency in Hz
T=1/f;%signal period
fs= 1000*f;% sampling frequency
ts= 1/fs;% sampling time
A = 0.2;% signal amplitude
%---------------------
sigma=[0 0.1:0.1:0.8 1:10];
K = 1001; %number of period
M = 9;% number of simulation times
TotalT=K*T;% simulation time
signal=A*sin(2*pi*f*(0:ts:TotalT))';
signal=signal+0.3*randn(size(signal));
ArrayN=1;
for p=1:1:length(sigma)

    fprintf('The noise index P=%d.n',p);
    Nsigma=sigma(p)% noise power density 2D=sigma^2 ts,sigma^2 is noise variance

    Ryy1tau=zeros(1001,M+1);

    Ryy2tau=zeros(1001,M+1);

    meany1t=zeros(1001,M+1);

    meany2t=zeros(1001,M+1);

%-------------------------------------------

  for m=0:M

      fprintf('Simulation time=%d.n',m);

      y1out=repmat(signal,1,ArrayN)+Nsigma*(2*rand(TotalT/ts+1,ArrayN)-1);

      y2out=sum(tanh(5*y1out),2);

   %-------------------E[y(t)]-----------

   %y1 = zeros(1001,1);

   y2 = zeros(1001,1);
   
   %y1=sum(reshape(y1out(1001:end-1),[1000,K-1]))/(K-1);
   %y2=sum(reshape(y2out(1001:end-1),[1000,K-1]))/(K-1);
   %y1=[y1;y1(1)];
   %y2=[y2';y2(1)];
   for k=1:1:K-1

      %y1=y1+y1out(k*1000+1:(k+1)*1000+1);

      y2=y2+y2out(k*1000+1:(k+1)*1000+1);

   end

   %y1=y1/(K-1);% E[y(t)] in [0 Ts[

   y2=y2/(K-1);

   %meany1t(:,m+1)=y1;
   meany2t(:,m+1)=y2;

%-----------------------------Ryy[tau]

    %y1out=y1out(1001:length(y1out));
    y2out=y2out(1001:length(y2out));

    %--------

    for n=0:1:1000

        tau=n*ts;    

        %y1tau=y1out(1:length(y1out)-n).*y1out(n+1:length(y1out));%y(t)y(t+tau)

        y2tau=y2out(1:length(y2out)-n).*y2out(n+1:length(y2out));%y(t)y(t+tau)

        %Eyy1=zeros(1001,1);

        Eyy2=zeros(1001,1);

          for k=0:1:K-3

              %Eyy1=Eyy1+y1tau(k*1000+1:(k+1)*1000+1);

              Eyy2=Eyy2+y2tau(k*1000+1:(k+1)*1000+1);

          end

        %Eyy1=Eyy1/(K-2);

        Eyy2=Eyy2/(K-2);      

        %Ryy1tau(n+1,m+1)=sum(Eyy1)/1000;

        Ryy2tau(n+1,m+1)=sum(Eyy2)/1000;

    end

end

%meany1 = sum(meany1t')/(M+1);

meany2 = sum(meany2t')/(M+1);

%Ryy1 = sum(Ryy1tau')/(M+1);

Ryy2 = sum(Ryy2tau')/(M+1);

%meany1=[meany1 meany1(2:1001)];

meany2=[meany2 meany2(2:1001)];

 for n=0:1000

   %EyEy1=meany1(1:length(meany1)-n).*meany1(n+1:length(meany1));%E(y(t))E(y(t+tau))

   %EyEy1=EyEy1(1:1001)';

   EyEy2=meany2(1:length(meany2)-n).*meany2(n+1:length(meany2));%E(y(t))E(y(t+tau))

   EyEy2=EyEy2(1:1001)';

   %Cyy1(n+1)=Ryy1(n+1)-sum(EyEy1)/1000;

   Cyy2(n+1)=Ryy2(n+1)-sum(EyEy2)/1000;

    %if Cyy2(n+1)<=0

     %  break

   %end

  end

 %Cyy1(n+2:1001)=0;

 Cyy2(n+2:1001)=0;

 %F1=(Cyy1(1)+2*sum(Cyy1(2:1001).*exp(-2*pi*i*(0.001:0.001:1))))*0.1;

 F2=(Cyy2(1)+2*sum(Cyy2(2:1001).*exp(-2*pi*i*(0.001:0.001:1))))*0.1;

 %Y1=sum(meany1(1:1001).*exp(-2*pi*i*(0:0.001:1)))/1000;

 Y2=sum(meany2(1:1001).*exp(-2*pi*i*(0:0.001:1)))/1000;

 SNRout(p)=abs(Y2)^2/real(F2)/f

%save 'SNRsinein' SNRout;

end