Theory and Application of Digital Signal Processing, Lawrence R.Rabiner-Bernard Gold, Prentice-Hall

メモ：いま勉強中です。

```r

#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})$である。

```r

#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})$である。

```r

#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)

```