function [At Ph A2]=EnvelopPhase(X,f0,fs)
HX=imag(hilbert(X));
[M N]=size(X);
t=0:1/fs:((N-1)/fs);
Ac=X.*cos(2*pi*f0*t)+HX.*sin(2*pi*f0*t);
As=HX.*cos(2*pi*f0*t)-X.*sin(2*pi*f0*t);
Ph=atan(As./Ac);
A2=Ac.*Ac+As.*As;
At=sqrt(A2);

N=10000;f0=10000;deltf=400;fs=22000;M=50;
a1=2;a2=4;a3=8; %（1）
sit1=pi/6;sit2=pi/4;sit3=pi/3;
X=Narrowbandsignal(N,f0,deltf,fs,M); %（2）
X=X/sqrt(var(X));
t=0:1/fs:((N-1)/fs);
X1=X+a1*cos(2*pi*f0*t+sit1);
X2=X+a2*cos(2*pi*f0*t+sit2);
X3=X+a3*cos(2*pi*f0*t+sit3);
[At1 Ph1 A21]=EnvelopPhase(X1,f0,fs);
[At2 Ph2 A22]=EnvelopPhase(X2,f0,fs);
[At3 Ph3 A23]=EnvelopPhase(X3,f0,fs);LA=0:0.4:12;
GA1=hist(At1,LA); %（3）
GA2=hist(At2,LA);
GA3=hist(At3,LA);
plot(LA,GA1,'--',LA,GA2,'-',LA,GA3,'-');
figure;
LP=-pi/2:0.05:pi/2;
GP1=hist((Ph1-sit1),LP);
GP2=hist((Ph2-sit2),LP);
GP3=hist((Ph3-sit3),LP);
plot(LP,GP1,'--',LP,GP2,'-',LP,GP3,'-');
figure;
LA2=0:1:120;
GA21=hist(A21,LA2);
GA22=hist(A22,LA2);
GA23=hist(A23,LA2);
plot(LA2,GA21,'--',LA2,GA22,'-',LA2,GA23,'-');
N=10000;f0=10000;deltf=400;fs=22000;M=50;