ftz = [-3.7431j -1.8051j 1.5699j 6.1910j];
RL = 22;
N = 4;
ftz = ftz/j;
nz = length(ftz);
syms w;
U = w-1/ftz(1);
V = ((w^2-1)^0.5)*(1-1/(ftz(1)^2))^0.5;
for k=2:1:N
PreU = U;
PreV = V;
if k>nz
U = CalU(inf, PreU, PreV);
V = CalV(inf, PreU, PreV);
else
U= CalU(ftz(k), PreU, PreV);
V = CalV(ftz(k), PreU, PreV);
end
end
F=sym2poly(U);
frz = roots(F);
P = poly(ftz);
F = poly(frz);
rip = 1./sqrt(10^(0.1*RL)-1.0)*abs(polyval(P,1)/polyval(F,1));
PP = conv(P,P);
FF = rip^2*conv(F,F);
EE = [zeros(1,length(FF)-length(PP)),PP]+FF;
r = roots(EE);
r = r(find(imag(r)>0));
E = poly(j*r);
F = poly(j*frz);
P = poly(j*ftz);
if mod(N-nz,2)==0
P = j*P;
end
EF = E+F;
m1 = zeros(1,N+1);
n1 = zeros(1,N+1);
for k=N+1:-2:1
n1(k) = j*imag(EF(k));
m1(k) = real(EF(k));
end
for k=N:-2:1
m1(k) = j*imag(EF(k));
n1(k) = real(EF(k));
end
if nz==N
epr = rip/sqrt(rip^2-1);
msl = epr/rip/(epr+1);
else
epr = 1.0;
msl = 0.0;
end
y21n = P/rip;
if mod(N,2)
if nz==N
y21n = y21n - j*msl*n1;
end
[r21,eigval,R] = residue(y21n,n1);
[r22,eigval,R] = residue(m1,n1);
else
if nz==N
y21n = y21n - j*msl*m1;
end
[r21,eigval,R] = residue(y21n,m1);
[r22,eigval,R] = residue(n1,m1);
end
r21 = real(r21);
r22 = real(r22);
Tnk = sqrt(r22);
T1k = r21./Tnk;
R1 = sum(T1k.^2);
RN = sum(Tnk.^2);
M = -diag([0;imag(eigval);0]);
M(1,2:N+1) = T1k.';
M(N+2,2:N+1) = Tnk.';
M(2:N+1,1) = T1k;
M(2:N+1,N+2) = Tnk;
M(1,N+2) = M(1,N+2)+msl;
M(N+2,1) = M(N+2,1)+msl;
M = round(M*10000)/10000
R1 = 1;
RN = 1;
w1 = -8;
w2 = 8;
dw = 0.01;
w = w1:dw:w2;
S21 = zeros(1, length(w));
S11 = zeros(1, length(w));
Tg = zeros(1, length(w));
for k=1:1:length(w)
R=zeros(N+2);
R(1,1)=R1;
R(N+2,N+2)=RN;
U = eye(N+2);
U(1, 1) = 0;
U(N+2, N+2) = 0;
Z = w(k)*U-j*R+M;
Zt = inv(Z);
S21(k) = 20*log10(-2*j*sqrt(R1*RN)*Zt(N+2,1));
S11(k) = 20*log10(1+2*j*R1*Zt(1,1));
for kk=2:N+1
Tg(k)=Tg(k)+(Zt(N+2,kk)*Zt(kk,1)/Zt(N+2,1));
end
Tg(k)=imag(Tg(k));
end
figure(1)
grid on
plot(w,real(S21),'g',w,real(S11),'b');
legend('S21','S11',2);
xlabel('归一化频率(Hz)');
ylabel('衰减(db)')
figure(2)
grid on
plot(w,Tg,'r');
legend('群时延特性');
xlabel('归一化频率(Hz)');
ylabel('群时延(ns)')