%Gaussian simulation of signal detection and width
Dv=0.1;%set spectral resolution
wc=-6.4:Dv:6.4;Np=128;%spectral velocity domain (128 pt fft
dbn=-30;Sn0=10.^(dbn/10);%specify noise level
Ax=7.5;A=10^((dbn+Ax)/10);%set signal strenth
sigc=.8;%set signal width
dbs=A/sqrt(2*pi)./sigc./exp((wc).^2./2./sigc^2);%compute gaussian signal distribution
figure;hold on;
for i=1:2000;
Sn=Sn0*(1+1/sqrt(8).*randn([1,length(wc)]));%add random noise for 8 spectral averages
Sp=dbs.*(1+1/sqrt(8).*randn([1,length(wc)]))+Sn;%add noise to signal
Sp2=dbs+Sn0+sqrt(dbs.^2+Sn0.^2)/sqrt(8).*randn([1,length(wc)]);
%figure;plot(wc,10*log10(Sp),'--',wc,10*log10(Sn));
plot(wc,10*log10(Sp)-10*log10(Sn),wc,10*log10(Sn)-10*log10(mean(Sn)));
%compute integrals
real([sum(Sp*.1) sum(Sn*.1) 10*log10(sum(Sp*.1)) 10*log10(sum(Sn*.1)) 10*log10(sum((Sp-Sn)*.1))  10*log10(sum((Sp-Sn)*.1)/sum(Sn*.1))] );
Spi(i,:)=Sp;Sni(i,:)=Sn;
%Find pts where signal exceed noise+sigma and noise + 2.5 sigma
ii=find(Sp>Sn*(1+1./sqrt(8)));jj=find(Sp>Sn*(1+2.5./sqrt(8)));
jji(i)=length(ii);
Snr= 10*log10(sum((Sp-Sn)*.1)/sum(Sn*.1));%snr is integral of signal/integral of noise
Snr0=10*log10(sqrt(2*pi)*2.5./Np/sqrt(8));
wbar=sum(wc.*(Sp-Sn)*.1)/sum((Sp-Sn)*.1);
sigmawidth=sqrt(sum((wc-wbar).^2.*(Sp-Sn)*.1)/sum((Sp-Sn)*.1));
%Compute V-Vmin
vv=(Snr-Snr0)/4.34+log(Dv*2.5/sigc);
Dvnor=sqrt(2*vv);
[Snr Dv*length(ii)/sigc/2 length(jj) Dvnor sigmawidth];

%Sp=Sp2;
Z=sum((Sp-Sn)*.1);Z0=sum((Sn)*.1);
pk=max(Sp(2:Np-1));pkn=max(Sn(2:Np-1));
if pk>=pkn
    ij=find(Sp>pk*.99);
    mx=find(ij>1 & ij<Np);
    ik=ij(mx(1));
    kk=ik;
    while (Sp(kk)>Sn0 );
       kk=kk-1;
       if kk==1;
           break;
       end;
    end;
    wn(i)=mean(wc(kk:kk+1));
    ki=ik;
    while (Sp(ki)>Sn0);
       ki=ki+1;
       if ki==Np;
           break;
       end;
    end;
    wx(i)=mean(wc(ki-1:ki));
    kz=kk+1:ki;
    Z2=sum((Sp(kz)-Sn0)*.1);
    Snr2= 10*log10(Z2/Z0);%snr is integral of signal/integral of noise
    wbar=sum(wc(kz).*(Sp(kz)-Sn0)*.1)/Z;
    width=(sum((wc(kz)-wbar).^2.*(Sp(kz)-Sn0)*.1)/Z);

else
    ij=find(Sn>pkn*.9);
    ik=ij(1);
    kk=ik;
    while Sn(kk)>Sn0
       kk=kk-1;
    end;
    wn(i)=mean(wc(kk:kk+1));
    ki=ik;
    while Sn(ki)>Sn0
       ki=ki+1;
    end;
    wx(i)=mean(wc(ki-1:ki));
    kz=kk+1:ki;
    Z2=sum((Sn(kz)-Sn0)*.1);
    Snr2= 10*log10(Z2/Z0);%snr is integral of signal/integral of noise
    wbar=sum(wc(kz).*(Sn(kz)-Sn0)*.1)/Z;
    width=(sum((wc(kz)-wbar).^2.*(Sn(kz)-Sn0)*.1)/Z);
end;

[Snr2 wbar length(jj) ki-kk sigmawidth];
swx(i)=width;
Zz(i)=Z2;
wb(i)=wbar;
kki(i)=ki-kk-1;
Snr2x(i)= Snr2;
if Zz(i)<0
      Zz(i)=Z;
end;
end;

[Ax Snr 10*log10(mean(Zz)/Z0) std(wb) sqrt(median((swx)))  std(sqrt(swx)) mean(kki) mean(jji)]
[Snr Dv*length(ii)/sigc/2 length(jj) Dvnor sigmawidth]
