| mehrfache Ausführung für verschieden starkes Rauschen |
N=64
dt=0.1
T0=2.8
tau=5.0
fs=2.6
phi=0.4
amp=10.0
t=arange(0,N)*dt
NN=20 # verschiedene Rauschleistungen
MM=1000 # Mittelungen
np=zeros(NN)
m1=zeros(NN)
m2=zeros(NN)
m3=zeros(NN)
s1=zeros(NN)
s2=zeros(NN)
s3=zeros(NN)
np=zeros(NN)
for i in range (0,NN):
np[i]=exp(i-10);
for j in range(0,MM):
x=amp*exp(-8*(t-T0)**2/tau**2)*cos(2*pi*fs*t+phi)+normal(0,sqrt(np[i]),N)
# 1. 3-Punkt-Schätzer
X=fft(append(x,zeros(N)))
df=1/float(2*N*dt)
E=abs(X)**2*dt**2
M=argmax(abs(E[0:N]))
Y=log(abs(E))
d=(Y[M-1]-Y[M+1])/(2*(Y[M-1]-2*Y[M]+Y[M+1]))
fe1=(M+d)*df
# 2. iterativer, komplexer Fit
fe2=fe1
for k in range(0,15):
C=cos(2*pi*fe2*t)
S=sin(2*pi*fe2*t)
R2=(sum(C*x)**2+sum(S*x)**2)/float(N**2)
dR2df=4*pi*(sum(S*x)*sum(t*C*x)-sum(C*x)*sum(t*S*x))/float(N**2)
d2R2df2=8*pi**2*(sum(S*x)**2-sum(C*x)*sum(t**2*C*x)+sum(C*x)**2-sum(S*x)*sum(t**2*S*x))/float(N**2)
fe2-=dR2df/d2R2df2
# 3. iterativer, komplexer Fit mit Gewichtung
hx=hilbert(x)
g=abs(hx)
fe3=fe1
for k in range(0,15):
C=cos(2*pi*fe3*t)
S=sin(2*pi*fe3*t)
R2=(sum(C*g*x)**2+sum(S*g*x)**2)/float(sum(g)**2)
dR2df=4*pi*(sum(S*g*x)*sum(t*C*g*x)-sum(C*g*x)*sum(t*S*g*x))/float(sum(g)**2)
d2R2df2=8*pi**2*(sum(S*g*x)**2-sum(C*g*x)*sum(t**2*C*g*x)+sum(C*g*x)**2-sum(S*g*x)*sum(t**2*S*g*x))/float(sum(g)**2)
fe3-=dR2df/d2R2df2
# Ergebnisse summieren
m1[i]+=fe1
m2[i]+=fe2
m3[i]+=fe3
s1[i]+=fe1**2
s2[i]+=fe2**2
s3[i]+=fe3**2
m1[i]/=MM
m2[i]/=MM
m3[i]/=MM
s1[i]=MM*(s1[i]/MM-m1[i]**2)/(MM-1)
s2[i]=MM*(s2[i]/MM-m2[i]**2)/(MM-1)
s3[i]=MM*(s3[i]/MM-m3[i]**2)/(MM-1)
print(m1[i],s1[i])
print(m2[i],s2[i])
print(m3[i],s3[i])
semilogx(np,m1,np,m2,np,m3)
show()
loglog(np,s1,np,s2,np,s3,np,64*np*dt/(pi**(5/2.0)*amp**2*tau**3))
show()
|
N=64;
dt=0.1;
T0=2.8;
tau=5.0;
fs=2.6;
phi=0.4;
amp=10;
t=(0:N-1)*dt;
NN=20; % verschiedene Rauschleistungen
MM=1000; % Mittelungen
np=zeros(1,NN);
m1=zeros(1,NN);
m2=zeros(1,NN);
m3=zeros(1,NN);
s1=zeros(1,NN);
s2=zeros(1,NN);
s3=zeros(1,NN);
np=zeros(1,NN);
for i=1:NN
np(i)=exp(i-11);
for j=1:MM
x=amp*exp(-8*(t-T0).^2/tau^2).*cos(2*pi*fs*t+phi)+normrnd(0,sqrt(np(i)),[1,N]);
% 1. 3-Punkt-Schätzer
X=fft([x,zeros(1,N)]);
df=1/(2*N*dt);
E=abs(X).^2*dt^2;
[mx,M]=max(abs(E(1:N)));
Y=log(abs(E));
d=(Y(mod(M-2,N)+1)-Y(mod(M,N)+1))/(2*(Y(mod(M-2,N)+1)-2*Y(M)+Y(mod(M,N)+1)));
fe1=(M-1+d)*df;
% 2. iterativer, komplexer Fit
fe2=fe1;
for k=1:15
C=cos(2*pi*fe2*t);
S=sin(2*pi*fe2*t);
R2=(sum(C.*x)^2+sum(S.*x)^2)/N;
dR2df=4*pi*(sum(S.*x)*sum(t.*C.*x)-sum(C.*x)*sum(t.*S.*x))/N;
d2R2df2=8*pi^2*(sum(S.*x)^2-sum(C.*x)*sum(t.^2.*C.*x)+sum(C.*x)^2-sum(S.*x)*sum(t.^2.*S.*x))/N;
fe2=fe2-dR2df/d2R2df2;
end
% 3. iterativer, komplexer Fit
X=fft(x);
X(1)=0;
for k=2:fix((length(x)+1)/2)
X(k)=-1j*X(k);
end
for k=fix((length(x)+3)/2):fix((length(x)+2)/2)
X(k)=0;
end
for k=fix((length(x)+4)/2):length(x)
X(k)=+1j*X(k);
end
hx=x+1j*ifft(X);
g=abs(hx);
fe3=fe1;
for k=1:15
C=cos(2*pi*fe3*t);
S=sin(2*pi*fe3*t);
R2=(sum(C.*g.*x)^2+sum(S.*g.*x)^2)/sum(g)^2;
dR2df=4*pi*(sum(S.*g.*x)*sum(t.*C.*g.*x)-sum(C.*g.*x)*sum(t.*S.*g.*x))/sum(g)^2;
d2R2df2=8*pi^2*(sum(S.*g.*x)^2-sum(C.*g.*x)*sum(t.^2.*C.*g.*x)+sum(C.*g.*x)^2-sum(S.*g.*x)*sum(t.^2.*S.*g.*x))/sum(g)^2;
fe3=fe3-dR2df/d2R2df2;
end
% Ergebnisse summieren
m1(i)=m1(i)+fe1;
m2(i)=m2(i)+fe2;
m3(i)=m3(i)+fe3;
s1(i)=s1(i)+fe1^2;
s2(i)=s2(i)+fe2^2;
s3(i)=s3(i)+fe3^2;
end
m1(i)=m1(i)/MM;
m2(i)=m2(i)/MM;
m3(i)=m3(i)/MM;
s1(i)=MM*(s1(i)/MM-m1(i)^2)/(MM-1);
s2(i)=MM*(s2(i)/MM-m2(i)^2)/(MM-1);
s3(i)=MM*(s3(i)/MM-m3(i)^2)/(MM-1);
[m1(i),s1(i)]
[m2(i),s2(i)]
[m3(i),s3(i)]
end
semilogx(np,m1,np,m2,np,m3)
waitforbuttonpress
loglog(np,s1,np,s2,np,s3,np,64*np*dt/(pi^(5/2)*amp^2*tau^3))
|
| Vergleich eines Least-Mean-Square-Schätzers mit einem Maximum-Likelihood-Schätzer eines linearen Trends für Poisson-verteilte Werte |
def det2(c11,c12,c21,c22):
return(c11*c22-c12*c21)
lfa=[0.0]
def logfactorial(i):
if i>=len(lfa):
for j in range(len(lfa),i+1):
lfa.append(log(j)+lfa[j-1])
return(lfa[i])
def loglikelihood(cv):
ls=0.0
er=False
for i in range(0,int(N[k])):
if cv[0,0]+cv[1,0]*t[i]>0:
ls+=x[i]*log(cv[0,0]+cv[1,0]*t[i])-(cv[0,0]+cv[1,0]*t[i])-logfactorial(x[i])
else:
er=True
if er:
ls=inf
return(ls)
lambda0=10.0
nu=5.0
M=10000
K=8
N=zeros(K)
S1cLMS=zeros(K)
S1dLMS=zeros(K)
S2cLMS=zeros(K)
S2dLMS=zeros(K)
S1cML=zeros(K)
S1dML=zeros(K)
S2cML=zeros(K)
S2dML=zeros(K)
CRLBc=zeros(K)
CRLBd=zeros(K)
for k in range(0,K):
N[k]=k+3
print('N=',N[k])
t=arange(N[k])+1
I=zeros((2,2))
for i in range(int(N[k])):
I[0,0]+=1.0/(lambda0+nu*t[i])
I[0,1]+=1.0/(lambda0+nu*t[i])*t[i]
I[1,0]+=1.0/(lambda0+nu*t[i])*t[i]
I[1,1]+=1.0/(lambda0+nu*t[i])*t[i]**2
Iinv=inv(I)
CRLBc[k]=Iinv[0,0]
CRLBd[k]=Iinv[1,1]
for i in range(M):
x=poisson(lambda0+nu*t)
#LMS
c=2.0*(2*N[k]+1)/float(N[k]*(N[k]-1))*sum(x)-6.0/float(N[k]*(N[k]-1))*sum(t*x)
d=12.0/float(N[k]*(N[k]**2-1))*sum(t*x)-6.0/float(N[k]*(N[k]-1))*sum(x)
S1cLMS[k]+=c
S1dLMS[k]+=d
S2cLMS[k]+=c**2
S2dLMS[k]+=d**2
#ML
better=True
while better:
cold=array([[c],[d]])
dldc=sum(x/(c+d*t))-float(N[k])
dldd=sum(t*x/(c+d*t))-float(N[k]*(N[k]+1))/2.0
dldcc=-sum(x/((c+d*t)**2))
dldcd=-sum((t*x)/((c+d*t)**2))
dlddd=-sum((t**2*x)/((c+d*t)**2))
lp=array([[dldc],[dldd]])
J=array([[dldcc,dldcd],[dldcd,dlddd]])
cnew=cold-dot(inv(J),lp)
if loglikelihood(cnew)>loglikelihood(cold):
c=cnew[0,0]
d=cnew[1,0]
better=True
else:
better=False
S1cML[k]+=c
S1dML[k]+=d
S2cML[k]+=c**2
S2dML[k]+=d**2
plot(N,lambda0*ones(K),'b',N,S1cLMS/float(M),'or',N,S1cML/float(M),'og',N,nu*ones(K),'b',N,S1dLMS/float(M),'or',N,S1dML/float(M),'og')
show()
plot(N,CRLBc,'b',N,M/float(M-1)*(S2cLMS/float(M)-(S1cLMS/float(M))**2),'or',N,M/float(M-1)*(S2cML/float(M)-(S1cML/float(M))**2),'og',\
N,CRLBd,'b',N,M/float(M-1)*(S2dLMS/float(M)-(S1dLMS/float(M))**2),'or',N,M/float(M-1)*(S2dML/float(M)-(S1dML/float(M))**2),'og')
show()
|
Matlab:
%nicht im Command Window
lambda0=10.0;
nu=5.0;
M=10000;
K=8;
N=zeros(1,K);
S1cLMS=zeros(1,K);
S1dLMS=zeros(1,K);
S2cLMS=zeros(1,K);
S2dLMS=zeros(1,K);
S1cML=zeros(1,K);
S1dML=zeros(1,K);
S2cML=zeros(1,K);
S2dML=zeros(1,K);
CRLBc=zeros(1,K);
CRLBd=zeros(1,K);
for k=1:K
N(k)=k+2;
N(k)
t=(1:N(k));
I=zeros(2,2);
for i=1:N(k)
I(1,1)=I(1,1)+1.0/(lambda0+nu*t(i));
I(1,2)=I(1,2)+1.0/(lambda0+nu*t(i))*t(i);
I(2,1)=I(2,1)+1.0/(lambda0+nu*t(i))*t(i);
I(2,2)=I(2,2)+1.0/(lambda0+nu*t(i))*t(i).^2;
end
Iinv=inv(I);
CRLBc(k)=Iinv(1,1);
CRLBd(k)=Iinv(2,2);
x=zeros(1,N(k));
for i=1:M
for j=1:N(k)
x(j)=poissrnd(lambda0+nu*t(j));
end
%LMS
c=2.0*(2*N(k)+1)/(N(k)*(N(k)-1))*sum(x)-6.0/(N(k)*(N(k)-1))*sum(t.*x);
d=12.0/(N(k)*(N(k)^2-1))*sum(t.*x)-6.0/(N(k)*(N(k)-1))*sum(x);
S1cLMS(k)=S1cLMS(k)+c;
S1dLMS(k)=S1dLMS(k)+d;
S2cLMS(k)=S2cLMS(k)+c^2;
S2dLMS(k)=S2dLMS(k)+d^2;
%ML
better=true;
while better
cold=[c;d];
dldc=sum(x./(c+d*t))-N(k);
dldd=sum(t.*x./(c+d*t))-(N(k)*(N(k)+1))/2.0;
dldcc=-sum(x./((c+d*t).^2));
dldcd=-sum((t.*x)./((c+d*t).^2));
dlddd=-sum((t.^2.*x)./((c+d*t).^2));
lp=[dldc;dldd];
J=[dldcc dldcd;dldcd dlddd];
cnew=cold-inv(J)*lp;
if loglikelihood(cnew,N(k),t,x)>loglikelihood(cold,N(k),t,x)
c=cnew(1);
d=cnew(2);
better=true;
else
better=false;
end
end
S1cML(k)=S1cML(k)+c;
S1dML(k)=S1dML(k)+d;
S2cML(k)=S2cML(k)+c^2;
S2dML(k)=S2dML(k)+d^2;
end
end
plot(N,lambda0*ones(1,K),'b',N,S1cLMS/M,'or',N,S1cML/M,'og',N,nu*ones(1,K),'b',N,S1dLMS/M,'or',N,S1dML/M,'og')
waitforbuttonpress
%
plot(N,CRLBc,'b',N,M/(M-1)*(S2cLMS/M-(S1cLMS/M).^2),'or',N,M/(M-1)*(S2cML/M-(S1cML/M).^2),'og',
N,CRLBd,'b',N,M/(M-1)*(S2dLMS/M-(S1dLMS/M).^2),'or',N,M/(M-1)*(S2dML/M-(S1dML/M).^2),'og')
function res=det2(c11,c12,c21,c22)
res=c11*c22-c12*c21
end
function res=logfactorial(i)
s=0;
for j=1:i
s=s+log(j);
end
res=s;
end
function res=loglikelihood(cv,N,t,x)
ls=0.0;
er=false;
for i=1:N
if cv(1)+cv(2)*t(i)>0
ls=ls+x(i)*log(cv(1)+cv(2)*t(i))-(cv(1)+cv(2)*t(i))-logfactorial(x(i));
else
er=true;
end
end
if er
ls=Inf;
end
res=ls;
end
Octave:
function res=det2(c11,c12,c21,c22)
res=c11*c22-c12*c21
end
function res=logfactorial(i)
s=0;
for j=1:i
s=s+log(j);
end
res=s;
end
function res=loglikelihood(cv,N,t,x)
ls=0.0;
er=false;
for i=1:N
if cv(1)+cv(2)*t(i)>0
ls=ls+x(i)*log(cv(1)+cv(2)*t(i))-(cv(1)+cv(2)*t(i))-logfactorial(x(i));
else
er=true;
end
end
if er
ls=Inf;
end
res=ls;
end
lambda0=10.0;
nu=5.0;
M=10000;
K=8;
N=zeros(1,K);
S1cLMS=zeros(1,K);
S1dLMS=zeros(1,K);
S2cLMS=zeros(1,K);
S2dLMS=zeros(1,K);
S1cML=zeros(1,K);
S1dML=zeros(1,K);
S2cML=zeros(1,K);
S2dML=zeros(1,K);
CRLBc=zeros(1,K);
CRLBd=zeros(1,K);
for k=1:K
N(k)=k+2;
N(k)
t=(1:N(k));
I=zeros(2,2);
for i=1:N(k)
I(1,1)=I(1,1)+1.0/(lambda0+nu*t(i));
I(1,2)=I(1,2)+1.0/(lambda0+nu*t(i))*t(i);
I(2,1)=I(2,1)+1.0/(lambda0+nu*t(i))*t(i);
I(2,2)=I(2,2)+1.0/(lambda0+nu*t(i))*t(i).^2;
end
Iinv=inv(I);
CRLBc(k)=Iinv(1,1);
CRLBd(k)=Iinv(2,2);
x=zeros(1,N(k));
for i=1:M
for j=1:N(k)
x(j)=poissrnd(lambda0+nu*t(j));
end
%LMS
c=2.0*(2*N(k)+1)/(N(k)*(N(k)-1))*sum(x)-6.0/(N(k)*(N(k)-1))*sum(t.*x);
d=12.0/(N(k)*(N(k)^2-1))*sum(t.*x)-6.0/(N(k)*(N(k)-1))*sum(x);
S1cLMS(k)=S1cLMS(k)+c;
S1dLMS(k)=S1dLMS(k)+d;
S2cLMS(k)=S2cLMS(k)+c^2;
S2dLMS(k)=S2dLMS(k)+d^2;
%ML
better=true;
while better
cold=[c;d];
dldc=sum(x./(c+d*t))-N(k);
dldd=sum(t.*x./(c+d*t))-(N(k)*(N(k)+1))/2.0;
dldcc=-sum(x./((c+d*t).^2));
dldcd=-sum((t.*x)./((c+d*t).^2));
dlddd=-sum((t.^2.*x)./((c+d*t).^2));
lp=[dldc;dldd];
J=[dldcc dldcd;dldcd dlddd];
cnew=cold-inv(J)*lp;
if loglikelihood(cnew,N(k),t,x)>loglikelihood(cold,N(k),t,x)
c=cnew(1);
d=cnew(2);
better=true;
else
better=false;
end
end
S1cML(k)=S1cML(k)+c;
S1dML(k)=S1dML(k)+d;
S2cML(k)=S2cML(k)+c^2;
S2dML(k)=S2dML(k)+d^2;
end
end
plot(N,lambda0*ones(1,K),'b',N,S1cLMS/M,'or',N,S1cML/M,'og',N,nu*ones(1,K),'b',N,S1dLMS/M,'or',N,S1dML/M,'og')
waitforbuttonpress
%
plot(N,CRLBc,'b',N,M/(M-1)*(S2cLMS/M-(S1cLMS/M).^2),'or',N,M/(M-1)*(S2cML/M-(S1cML/M).^2),'og',
N,CRLBd,'b',N,M/(M-1)*(S2dLMS/M-(S1dLMS/M).^2),'or',N,M/(M-1)*(S2dML/M-(S1dML/M).^2),'og')
|