非線形デジィタル信号処理　棟安・田口先生　朝倉書店

```r

#p.26 z変換の割り算
#p.84の例 書籍：z変換とその応用　柳井浩先生

h=100*c(1,0.8,0.6,0.6)
b=c(0.8,0.5,0.3)
K=length(h)-1
N=length(b)
H=array(0,dim=c(K+1,N+1))

for(j in 1:(K+1)){
hvec=h[j:1]
H[j,1:j]=hvec
}

a=H%*%c(1,b)

```

```r

#p.38 メジアンフィルタ

n=100
x=c(rnorm(n/2,-5,10),rnorm(n/2,5,10))
plot(x)

loglik_x=function(b,sig){
k=((sig^2)^(-b/2))*(gamma(3/b)*gamma(1/b))^(b/2)
return(-(n*log(b)+n*log(k)/b-sum(k*abs(x)^b)-n*log(2*gamma(1/b))))
}

#最尤法

beta=1
sigma=1
eta=10^(-5)
h=10^(-5)
ite=10^5


for(l in 1:ite){
beta=beta-eta*(loglik_x(beta+h,sigma)-loglik_x(beta,sigma))/h
sigma=sigma-eta*(loglik_x(beta,sigma+h)-loglik_x(beta,sigma))/h
#print(loglik_x(beta,sigma))
}

#推定されたパラメータの一般化ガウス分布
gauss=function(y){
k=(sigma^(-beta))*(gamma(3/beta)*gamma(1/beta))^(beta/2)
return((beta*k^(1/beta))*exp(-k*abs(y)^beta)/(2*gamma(1/beta)))
}

#MED密度関数
loglik_med_density=function(x,M){
dis=integrate(gauss,-Inf,x)$value
value=sum(log(c(1:n)))-2*n*sum(log(gamma(M+1)))+sum(log(gauss(x)))+M*sum(log(dis)+log(1-dis))
return(value)    
}

#MED分布関数
med_dis=function(x,M){
dis=integrate(gauss,-Inf,x)$value
return(1-pbinom(M,n,dis))
}

#MED密度関数・MED分布関数

#密度関数出力(対数尤度)
m=10
for(s in 1:m){
log_density_med=c()
for(j in 1:n){log_density_med=c(log_density_med,(loglik_med_density(x[j],s)))}
plot(x,log_density_med,main=paste0("M=",s))
}

#分布関数出力
J=5
m=n-3
for(j in 1:J){
dis_med_vec=c()
for(s in 1:m){dis_med_vec=c(dis_med_vec,med_dis(x[j],s))}
plot(c(1:m),dis_med_vec,main=paste0("x=",x[j]))
}

```

```r

#p.108 傾斜適応系εフィルタ

data("AirPassengers")
x=c(AirPassengers)
K=3
L=2
xchil=x[(L+1):(length(x)-K)]

#目的変数：入力x(n)
Y=c()
for(j in 1:length(xchil)){
Y=c(Y,rep(xchil[j],L+K+1))
}

#予測値x(n-k)+a*k, k=-K~L
f=function(a){
const=a[1]*c(-K:L)
pred=c()
for(j in (L+1):(L+length(xchil))){
pred=c(pred,(x[(j-L):(j+K)]+const))  
}
return(pred)
}

#ガウスニュートン法(wは全て1)

h=0.001
lam=1
ite=10^3
#sita=rep(0.01,1+L+K+1)
sita=0.01
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}
pred=f(sita)
plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),main="",col=2)
par(new=T)
plot(c(1:length(Y)),pred,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),main="Y:Pred",col=3)


```


```r

#2次フィルタ
rt=c(62.7,75.7,92.2,98.3,106.2,109.5,114.8)
Qt=rt*1.5+rnorm(length(rt))*rt

t=c(1:length(rt))

M=3

f=function(x){

h2=matrix(x,nrow=M)
Qhat=c()
for(j in 1:(length(rt)-M+1)){
  rtsub=rt[j:(M+j-1)]
  Qhat=c(Qhat,as.numeric(c(rtsub)%*%h2%*%c(rtsub)))
  }
return(Qhat[-length(Qhat)])
}


Y=Qt[(M+1):length(t)]
h=0.001
lam=1
ite=10^4
sita=rep(0.01,M^2)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}





```


```r

#p.128 双線形モデル
xt=c(62.7,75.7,92.2,98.3,106.2,109.5,114.8,118.9,120)
yt=sin(xt)+rnorm(length(xt))*xt

t=c(1:length(xt))

M=3

f=function(x){
h0=x[1:(M-1)];x=x[-c(1:(M-1))]
h1=x[1:M];x=x[-c(1:M)]
h2=matrix(x,nrow=M)
Qhat=c()
for(j in 1:(length(xt)-M+1)){
  xtsub=xt[j:(M+j-1)]
  ytsub=yt[j:(M+j-2)]
  Qhat=c(Qhat,sum(xtsub*h1)+sum(ytsub*h0)+as.numeric(c(xtsub)%*%h2%*%c(ytsub)))
  }
return(Qhat[-length(Qhat)])
}


Y=yt[(M+1):length(t)]
h=0.001
lam=1
ite=10^3
sita=rep(0.01,M+(M-1)+M*(M-1))
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}



```


```r

#LMS2次フィルタ
xt=c(62.7,75.7,92.2,98.3,106.2,109.5,114.8)
dn=c()
for(j in 1:length(xt)){
dn=c(dn,xt[j]*xt[j:(3+j-1)]*2+rnorm(1))  
}

t=c(1:length(xt))

M=3

f=function(x){
dhat=c()
for(j in 1:(length(xt)-M+1)){
  x2=x
  xtsub=xt[j:(M+j-1)]
  pred=c()
  for(i in 1:length(xtsub)){
    hsub=x2[1:((M-i+1))];x2=x2[-c(1:((M-i+1)))]
    pred=c(pred,sum(hsub*xtsub[1:(M-i+1)])*xtsub[(M-i+1)])
  }
  dhat=c(dhat,as.numeric(sum(pred)))
  }
return(dhat)
}


Y=dn[(M):length(t)]
h=0.001
lam=1
ite=10^4
sita=rep(0.01,M*(M+1)/2)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}





```