デジィタルフィルタ ー原理と設計法ー 科学情報出版 陶山健仁先生

r=function(n,a){
return(ifelse(n>=0,1,0)*a^n)
}

N=c(0:10)

#指数信号p.10
plot(N,r(N,0.5)) #a>0,a<1
plot(N,r(N,-0.5)) #a>0,a<1
plot(N,r(N,2)) #a>1
plot(N,r(N,-2)) #a<-1

#p.21

fs=1000
Ts=1/fs
x=rep(1,4)
n=c(1:length(x))-1
W=seq(-fs*10*0,fs*10,1)
w0=2*pi/(length(x)*Ts)

f=function(w){
mat=matrix(rep(w,length(n)),nrow=length(w))
values=apply(exp(-t(t(mat)*n*Ts*1i)),1,sum)
return(values)
}

Xw=sqrt(Re(f(W))^2+Im(f(W))^2)

#振幅スペクトル|X(w)|
plot(W,Xw)

#基本周波数の確認
plot(W[round(w0):(round(w0)*3)],Xw[round(w0):(round(w0)*3)])

#p.23

fs=1
Ts=1/fs
x=c(2,1,2,2,-1,-1,1,-2)
n=c(1:length(x))-1
W=seq(-20*fs*0/2,7,1)
w0=2*pi/(length(x)*Ts)

f=function(w){
mat=matrix(rep(w,length(n)),nrow=length(w))
values=apply(exp(-(t(mat)*(2*pi*n)*Ts*1i/length(x)))*x,2,sum)
return(values)
}

Xw=sqrt(Re(f(W))^2+Im(f(W))^2)

#振幅スペクトル|X(w)| 図1-2-19
plot(W,sqrt(Re(f(W))^2+Im(f(W))^2),type="h")

#p.37 インパルス応答
#複数の場合
dat=EuStockMarkets

#y:フィードバック成分として見る
y=dat[,1]

M=3 #自己回帰
N=4 #インパルス応答
S=max(M,N)

Y=y[1:length(y)]
X=array(0,dim=c(length(Y),M+N))

for(i in 1:(nrow(X))){
x=rep(0,N)
if(i<=N){x[i]=1} 

if(i<=(M+1)){

if(i==1){
X[i,]=c(rep(0,M),x)  #x(0)=delta(0)=1  
}else{
X[i,]=c(rep(0,M-i+1),-y[1:(i-1)],x)  #x(0)=delta(0)=1
}

}else{
X[i,]=c(-y[(i-M):(i-1)],x) #x(n)=delta(n)=0,n>0
}

}

#最小二乗法による推定
param=solve(t(X)%*%X)%*%t(X)%*%Y
bk=param
hn=X%*%param

#精度
cor(Y,X%*%param)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),col=2)
par(new=T)
plot(c(1:length(Y)),X%*%param,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),col=3)

#p.36 Ⅱ-1-3  線形時不変離散時間システム
dat=EuStockMarkets

#y:フィードバック成分として見る
y=dat[,1]

#x:入力信号としてみる
x=dat[,2]

M=3 #自己回帰
N=4 #入力信号部分
S=max(M,N)

Y=y[(S+1):length(y)]
X=array(0,dim=c(length(Y),M+N))

for(i in 1:nrow(X)){
X[i,]=c(-y[i:(M+i-1)],x[i:(N+i-1)]) 
}

#最小二乗法による推定
param=solve(t(X)%*%X)%*%t(X)%*%Y
bk=param[1:M]
ak=param[-c(1:M)]

#精度
cor(Y,X%*%param)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),col=2)
par(new=T)
plot(c(1:length(Y)),X%*%param,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),col=3)

#p.176 直線位相FIRフィルタ
data("LakeHuron")
w=c(1:length(c(LakeHuron)))
t=mean(w/(2*pi))
w0=2*pi/(length(w)*t)
w=w*w0

Hd=c(LakeHuron)

#ガウスニュートン法

M=5
X=t(t(w))%*%t(c(0:M))
f=function(x){
#重み  
return((cos(X)%*%x))  
}



Y=(Hd)
h=0.001
lam=1
ite=10^4

sita=solve(t(cos(X))%*%cos(X)+diag(lam,M+1))%*%t(cos(X))%*%Y
eta=10^(-1)

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
}

Wy=exp(X%*%sita)
sita=sita+eta*solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))

print(sum((Y-f(sita))^2))　
}

plot(w,f(sita),xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=2)
par(new=T)
plot(w,Y,xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=3,main=paste0("重相関:",cor(Y,f(sita))))

#p.196 FIRフィルタの複素近似
data("LakeHuron")
w=c(1:length(c(LakeHuron)))
t=mean(w/(2*pi))
w0=2*pi/(length(w)*t)
w=w*w0

Hd=c(LakeHuron)

#ガウスニュートン法

N=30
X=t(t(w))%*%t(c(0:N))
f=function(x){
y=x[-c(1:(N+1))]  
x=x[1:(N+1)]
val=cos(X)%*%x-1i*sin(X)%*%y
return(sqrt(Re(val)^2+Im(val)^2))  
}



Y=(Hd)
h=0.001
lam=0.001
ite=5000

sita=rep(1,2*(N+1))
eta=10^(-0)

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
}

sita=sita+eta*solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
#print(sum((Y-f(sita))^2))　
}

plot(w,f(sita),xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=2)
par(new=T)
plot(w,Y,xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=3,main=paste0("重相関:",cor(Y,f(sita))))

#p.214 等リプル設計と交番定理
data("LakeHuron")
Hd=c(LakeHuron)
Y=Hd

#周期を最急降下法で求める
f=function(w0){
w=w0*c(1:length(Hd))  
M=length(w)-2
A=t(t(w))%*%t(c(0:M))
A=cbind(cos(A),(-1)^c(0:(M+1)))
sita=solve(A)%*%Hd
return(abs(sita[length(sita)]))}

WW0=5
h=0.001
ite=10^3
eta=10^(-7)
for(l in 1:ite){
WW0_pre=WW0  
df=-(f(WW0+h)-f(WW0))/h
WW0=WW0+eta*df 
print(f(WW0))
}

WW=WW0*c(1:length(Hd))  
M=length(WW)-2
AA=t(t(WW))%*%t(c(0:M))
AA=cbind(cos(AA),(-1)^c(0:(M+1)))
param=solve(AA)%*%Hd