メモ：いま勉強中です。
#p.136
#3.30 Remez Exchange Algorithm Design of Optimal FIR FIlters
w=c(-4.1 ,
-4.099375 ,
-4.09875 ,
-4.098125 ,
-4.0975 ,
-4.096875 ,
-4.09625 ,
-4.095625 ,
-4.095 ,
-4.094375 ,
-4.09375 ,
-4.093125 ,
-4.0925 ,
-4.091875 ,
-4.09125 ,
-4.090625 ,
-4.09 ,
-4.089375 ,
-4.08875 ,
-4.088125 ,
-4.0875 ,
-4.086875 ,
-4.08625 ,
-4.085625 ,
-4.085 ,
-4.084375 ,
-4.08375 ,
-4.083125 ,
-4.0825 ,
-4.081875 ,
-4.08125 ,
-4.080625 ,
-4.08 ,
-4.079375 ,
-4.07875 ,
-4.078125 ,
-4.0775 ,
-4.076875 ,
-4.07625 ,
-4.075625 ,
-4.075 ,
-4.074375 ,
-4.07375 ,
-4.073125 ,
-4.0725 ,
-4.071875 ,
-4.07125 ,
-4.070625 ,
-4.07 ,
-4.069375 ,
-4.06875 ,
-4.068125 ,
-4.0675 ,
-4.066875 ,
-4.06625 ,
-4.065625 ,
-4.065 ,
-4.064375 ,
-4.06375 ,
-4.063125 ,
-4.0625 ,
-4.061875 ,
-4.06125 ,
-4.060625 ,
-4.06 ,
-4.059375 ,
-4.05875 ,
-4.058125 ,
-4.0575 ,
-4.056875 ,
-4.05625 ,
-4.055625 ,
-4.055 ,
-4.054375 ,
-4.05375 ,
-4.053125 ,
-4.0525 ,
-4.051875 ,
-4.05125 ,
-4.050625 ,
-4.05 ,
-4.049375 ,
-4.04875 ,
-4.048125 ,
-4.0475 ,
-4.046875 ,
-4.04625 ,
-4.045625 ,
-4.045 ,
-4.044375 ,
-4.04375 ,
-4.043125 ,
-4.0425 ,
-4.041875 ,
-4.04125 ,
-4.040625 ,
-4.04 ,
-4.039375 ,
-4.03875 ,
-4.038125
)
w=max(w)-w
w=2*pi*w*10
Y=c(0.015997742 ,
0.015997742 ,
0.016329822 ,
0.015333584 ,
0.06580965 ,
-0.012893164 ,
-0.023519705 ,
0.057839745 ,
0.002050408 ,
0.012344869 ,
0.058835983 ,
0.002382487 ,
0.016329822 ,
0.00105417 ,
0.010684472 ,
0.023967648 ,
0.042564093 ,
-0.003262862 ,
-0.051414373 ,
0.038247061 ,
0.020314774 ,
0.036918744 ,
0.088723128 ,
0.016329822 ,
0.121266908 ,
0.056843507 ,
-0.139415411 ,
0.09337224 ,
-0.077648645 ,
-0.115505695 ,
0.068466286 ,
0.046216967 ,
0.041235776 ,
-0.036470801 ,
0.084406096 ,
0.088723128 ,
-0.033814166 ,
-0.003927021 ,
0.030609235 ,
-0.077980724 ,
0.013341107 ,
0.080089064 ,
-0.012561085 ,
-0.002930783 ,
0.021643092 ,
0.012344869 ,
0.048209443 ,
0.044888649 ,
-0.040787833 ,
0.043560331 ,
0.036918744 ,
-0.017542276 ,
0.039243299 ,
-0.021859308 ,
-0.016546038 ,
0.020978933 ,
0.014337345 ,
0.023303489 ,
-0.010236529 ,
0.047545284 ,
0.091047684 ,
-0.051082294 ,
0.012676948 ,
0.10399878 ,
-0.021859308 ,
-0.012229006 ,
0.014669425 ,
-0.077648645 ,
0.052858554 ,
0.21225666 ,
0.048541522 ,
-0.184578208 ,
0.068134206 ,
0.052858554 ,
-0.10753579 ,
-0.024183864 ,
0.008691996 ,
0.010020313 ,
0.027288441 ,
0.018654377 ,
0.041899935 ,
0.014669425 ,
0.002714567 ,
0.02263933 ,
0.038247061 ,
0.050866078 ,
0.01168071 ,
0.007363678 ,
-0.013889403 ,
0.003046646 ,
0.061492618 ,
-0.039791595 ,
0.010020313 ,
0.107651653 ,
-0.013225244 ,
0.015665663 ,
0.030277156 ,
-0.057059723 ,
0.050866078 ,
-0.006583656
)
D=Y
r=length(w)
K=3
W=function(ws){
return(1/(1+exp(-K*abs(cos(ws)))))
}
Weji=1/W(w)
mat=array(0,dim=c(r,r-1))
for(j in 1:(r-1)){
mat[,j]=cos(j*w)
}
mat=cbind(rep(1,r),mat,Weji)
#ガウス・ザイデル法 p.248
A=t(mat)%*%mat
C=t(mat)%*%D
R=C/diag(A)
b=A/c(diag(A))
diag(b)=0
mat=cbind(R,-b)
ite=10^3
x=rep(1,r+1)/10
eta=0.01
for(l in 1:ite){
x_pre=x
z=c(1,x)
x=x*(1-eta)+eta*mat%*%z
df=x-x_pre
print(sqrt(sum(Re(df)^2+Im(df)^2)))
}
alpha=x[1:r]
delta=x[(r+1)]
pre=c()
for(j in 1:r){
pre=c(pre,sum(alpha*cos(w[j]*c(0:(r-1)))))
}
delta
plot(w,pre,xlim=c(min(w),max(w)),ylim=c(min(c(pre,Y)),max(c(pre,Y))),xlab="w",ylab="H",col=2)
par(new=T)
plot(w,Y,xlim=c(min(w),max(w)),ylim=c(min(c(pre,Y)),max(c(pre,Y))),xlab="w",ylab="H",col=3)
$IIR$フィルターの$z$変換
$$ H(z)=A\prod_{i = 1}^{K}\frac{1+a_i z^{-1}+b_{i} z^{-2}}{1+c_i z^{-1}+d_{i} z^{-2}}$$
をデータ$H_d(e^{jw_{i}})$について二乗誤差(ただし、$j=\sqrt{-1}$)：
$$ Q(\theta)=\sum_{i=1}^{M}(|H(e^{jw_{i}})|-|H_d(e^{jw_{i}})|)^{2} $$
を最小にするようにガウス法にてパラメータ推定を行ってみた。
$\theta$はパラメータで$\theta=(A,a_{1},...,a_{K},b_{1},...,b_{K},c_{1},...,c_{K},d_{1},...,d_{K})$である。
#p.268 1.Minimum Mean Squared Error Design
#IIR Filter
w=c(-4.1 ,
-4.099375 ,
-4.09875 ,
-4.098125 ,
-4.0975 ,
-4.096875 ,
-4.09625 ,
-4.095625 ,
-4.095 ,
-4.094375 ,
-4.09375 ,
-4.093125 ,
-4.0925 ,
-4.091875 ,
-4.09125 ,
-4.090625 ,
-4.09 ,
-4.089375 ,
-4.08875 ,
-4.088125 ,
-4.0875 ,
-4.086875 ,
-4.08625 ,
-4.085625 ,
-4.085 ,
-4.084375 ,
-4.08375 ,
-4.083125 ,
-4.0825 ,
-4.081875 ,
-4.08125 ,
-4.080625 ,
-4.08 ,
-4.079375 ,
-4.07875 ,
-4.078125 ,
-4.0775 ,
-4.076875 ,
-4.07625 ,
-4.075625 ,
-4.075 ,
-4.074375 ,
-4.07375 ,
-4.073125 ,
-4.0725 ,
-4.071875 ,
-4.07125 ,
-4.070625 ,
-4.07 ,
-4.069375 ,
-4.06875 ,
-4.068125 ,
-4.0675 ,
-4.066875 ,
-4.06625 ,
-4.065625 ,
-4.065 ,
-4.064375 ,
-4.06375 ,
-4.063125 ,
-4.0625 ,
-4.061875 ,
-4.06125 ,
-4.060625 ,
-4.06 ,
-4.059375 ,
-4.05875 ,
-4.058125 ,
-4.0575 ,
-4.056875 ,
-4.05625 ,
-4.055625 ,
-4.055 ,
-4.054375 ,
-4.05375 ,
-4.053125 ,
-4.0525 ,
-4.051875 ,
-4.05125 ,
-4.050625 ,
-4.05 ,
-4.049375 ,
-4.04875 ,
-4.048125 ,
-4.0475 ,
-4.046875 ,
-4.04625 ,
-4.045625 ,
-4.045 ,
-4.044375 ,
-4.04375 ,
-4.043125 ,
-4.0425 ,
-4.041875 ,
-4.04125 ,
-4.040625 ,
-4.04 ,
-4.039375 ,
-4.03875 ,
-4.038125
)
w=max(w)-w
w=2*pi*w*10
Hd=c(0.015997742 ,
0.015997742 ,
0.016329822 ,
0.015333584 ,
0.06580965 ,
-0.012893164 ,
-0.023519705 ,
0.057839745 ,
0.002050408 ,
0.012344869 ,
0.058835983 ,
0.002382487 ,
0.016329822 ,
0.00105417 ,
0.010684472 ,
0.023967648 ,
0.042564093 ,
-0.003262862 ,
-0.051414373 ,
0.038247061 ,
0.020314774 ,
0.036918744 ,
0.088723128 ,
0.016329822 ,
0.121266908 ,
0.056843507 ,
-0.139415411 ,
0.09337224 ,
-0.077648645 ,
-0.115505695 ,
0.068466286 ,
0.046216967 ,
0.041235776 ,
-0.036470801 ,
0.084406096 ,
0.088723128 ,
-0.033814166 ,
-0.003927021 ,
0.030609235 ,
-0.077980724 ,
0.013341107 ,
0.080089064 ,
-0.012561085 ,
-0.002930783 ,
0.021643092 ,
0.012344869 ,
0.048209443 ,
0.044888649 ,
-0.040787833 ,
0.043560331 ,
0.036918744 ,
-0.017542276 ,
0.039243299 ,
-0.021859308 ,
-0.016546038 ,
0.020978933 ,
0.014337345 ,
0.023303489 ,
-0.010236529 ,
0.047545284 ,
0.091047684 ,
-0.051082294 ,
0.012676948 ,
0.10399878 ,
-0.021859308 ,
-0.012229006 ,
0.014669425 ,
-0.077648645 ,
0.052858554 ,
0.21225666 ,
0.048541522 ,
-0.184578208 ,
0.068134206 ,
0.052858554 ,
-0.10753579 ,
-0.024183864 ,
0.008691996 ,
0.010020313 ,
0.027288441 ,
0.018654377 ,
0.041899935 ,
0.014669425 ,
0.002714567 ,
0.02263933 ,
0.038247061 ,
0.050866078 ,
0.01168071 ,
0.007363678 ,
-0.013889403 ,
0.003046646 ,
0.061492618 ,
-0.039791595 ,
0.010020313 ,
0.107651653 ,
-0.013225244 ,
0.015665663 ,
0.030277156 ,
-0.057059723 ,
0.050866078 ,
-0.006583656
)
M=length(w)
K=10
a=rep(1,K)
b=rep(1,K)
c=rep(1,K)
d=rep(1,K)
#Gauss法
f=function(x){
A=x[1];x=x[-1]
ak=x[1:K];x=x[-c(1:K)]
bk=x[1:K];x=x[-c(1:K)]
ck=x[1:K];x=x[-c(1:K)]
dk=x
val=rep(1,M)
for(j in 1:K){
val=val*(1+ak[j]*exp(1i*w)^(-1)+bk[j]*exp(1i*w)^(-2))/(1+ck[j]*exp(1i*w)^(-1)+dk[j]*exp(1i*w)^(-2))
}
val=sqrt(Re(A*val)^2+Im(A*val)^2)
return(val)
}
Y=abs(Hd)
h=0.001
lam=1
ite=10^3
sita=c(1,a,b,c,d)
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))
}
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)
同じく、また別の$z$変換
$$ H(z)=A\prod_{i = 1}^{K}\frac{\sum_{l=0}^{P}a_{l}z^{-l}}{\sum_{l=0}^{P}b_{l}z^{-l}} \space \space \space ただし、a_{0}=1、b_{0}=1$$
をデータ$H_d(e^{jw_{i}})$について二乗誤差：
$$ Q(\theta)=\sum_{i=1}^{M}(|H(e^{jw_{i}})|-|H_d(e^{jw_{i}})|)^{2} $$
を最小にするようにガウス法にてパラメータ推定を行ってみた。
$\theta$はパラメータで$\theta=(A,a_{1},...,a_{P},b_{1},...,b_{P})$である。
#p.268
#Other z'transform
w=c(-4.1 ,
-4.099375 ,
-4.09875 ,
-4.098125 ,
-4.0975 ,
-4.096875 ,
-4.09625 ,
-4.095625 ,
-4.095 ,
-4.094375 ,
-4.09375 ,
-4.093125 ,
-4.0925 ,
-4.091875 ,
-4.09125 ,
-4.090625 ,
-4.09 ,
-4.089375 ,
-4.08875 ,
-4.088125 ,
-4.0875 ,
-4.086875 ,
-4.08625 ,
-4.085625 ,
-4.085 ,
-4.084375 ,
-4.08375 ,
-4.083125 ,
-4.0825 ,
-4.081875 ,
-4.08125 ,
-4.080625 ,
-4.08 ,
-4.079375 ,
-4.07875 ,
-4.078125 ,
-4.0775 ,
-4.076875 ,
-4.07625 ,
-4.075625 ,
-4.075 ,
-4.074375 ,
-4.07375 ,
-4.073125 ,
-4.0725 ,
-4.071875 ,
-4.07125 ,
-4.070625 ,
-4.07 ,
-4.069375 ,
-4.06875 ,
-4.068125 ,
-4.0675 ,
-4.066875 ,
-4.06625 ,
-4.065625 ,
-4.065 ,
-4.064375 ,
-4.06375 ,
-4.063125 ,
-4.0625 ,
-4.061875 ,
-4.06125 ,
-4.060625 ,
-4.06 ,
-4.059375 ,
-4.05875 ,
-4.058125 ,
-4.0575 ,
-4.056875 ,
-4.05625 ,
-4.055625 ,
-4.055 ,
-4.054375 ,
-4.05375 ,
-4.053125 ,
-4.0525 ,
-4.051875 ,
-4.05125 ,
-4.050625 ,
-4.05 ,
-4.049375 ,
-4.04875 ,
-4.048125 ,
-4.0475 ,
-4.046875 ,
-4.04625 ,
-4.045625 ,
-4.045 ,
-4.044375 ,
-4.04375 ,
-4.043125 ,
-4.0425 ,
-4.041875 ,
-4.04125 ,
-4.040625 ,
-4.04 ,
-4.039375 ,
-4.03875 ,
-4.038125
)
w=max(w)-w
w=2*pi*w*10
Hd=c(0.015997742 ,
0.015997742 ,
0.016329822 ,
0.015333584 ,
0.06580965 ,
-0.012893164 ,
-0.023519705 ,
0.057839745 ,
0.002050408 ,
0.012344869 ,
0.058835983 ,
0.002382487 ,
0.016329822 ,
0.00105417 ,
0.010684472 ,
0.023967648 ,
0.042564093 ,
-0.003262862 ,
-0.051414373 ,
0.038247061 ,
0.020314774 ,
0.036918744 ,
0.088723128 ,
0.016329822 ,
0.121266908 ,
0.056843507 ,
-0.139415411 ,
0.09337224 ,
-0.077648645 ,
-0.115505695 ,
0.068466286 ,
0.046216967 ,
0.041235776 ,
-0.036470801 ,
0.084406096 ,
0.088723128 ,
-0.033814166 ,
-0.003927021 ,
0.030609235 ,
-0.077980724 ,
0.013341107 ,
0.080089064 ,
-0.012561085 ,
-0.002930783 ,
0.021643092 ,
0.012344869 ,
0.048209443 ,
0.044888649 ,
-0.040787833 ,
0.043560331 ,
0.036918744 ,
-0.017542276 ,
0.039243299 ,
-0.021859308 ,
-0.016546038 ,
0.020978933 ,
0.014337345 ,
0.023303489 ,
-0.010236529 ,
0.047545284 ,
0.091047684 ,
-0.051082294 ,
0.012676948 ,
0.10399878 ,
-0.021859308 ,
-0.012229006 ,
0.014669425 ,
-0.077648645 ,
0.052858554 ,
0.21225666 ,
0.048541522 ,
-0.184578208 ,
0.068134206 ,
0.052858554 ,
-0.10753579 ,
-0.024183864 ,
0.008691996 ,
0.010020313 ,
0.027288441 ,
0.018654377 ,
0.041899935 ,
0.014669425 ,
0.002714567 ,
0.02263933 ,
0.038247061 ,
0.050866078 ,
0.01168071 ,
0.007363678 ,
-0.013889403 ,
0.003046646 ,
0.061492618 ,
-0.039791595 ,
0.010020313 ,
0.107651653 ,
-0.013225244 ,
0.015665663 ,
0.030277156 ,
-0.057059723 ,
0.050866078 ,
-0.006583656
)
M=length(w)
L=20
#Gauss法
f=function(x){
A=x[1];x=x[-1]
ak=x[1:(L)];x=x[-c(1:(L))]
bk=x
val1=1
val2=1
for(j in 1:L){
val1=val1+ak[j]*exp(w*1i)^(-j)
val2=val2+bk[j]*exp(w*1i)^(-j)
}
val=A*val1/val2
val=sqrt(Re(A*val)^2+Im(A*val)^2)
return(val)
}
Y=abs(Hd)
h=0.001
lam=1
ite=10^3
sita=c(1,rep(1,2*L))
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))
}
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)