function [P,ryy,rxy,rle,v,d,w,b,mmse,SNRmfb_in_dB,c]=teq(p,L,mu,delta,sigma,Ex_bar,filtertype) % programmed by Ryoulhee Kwak % % filtertype 1=FIR , 2= IIR(just for one pole filter) % p=channel impulse response % in case of FIR, p= 1+0.5D+2.3D^2->p=[1 0.5 2.3] % in case of IIR, p=1/(2.4-0.5D)-> p=[2.4 -0.5] % sigma =>noise power in linear scale % delta =>delay % L => number of taps % c=> w*p % v=>eigen vectors % d=>eigent values %mu => for b if filtertype==1% FIR norm_p=norm(p); [temp size_p]=size(p); P=toeplitz([p(1), zeros(1,L-1)]',[p,zeros(1,L-1)]); ryy=Ex_bar*P*P'+sigma*eye(L,L); rxy=[zeros(mu+1,delta) eye(mu+1) zeros(mu+1,L+size_p-2-mu-delta)]*P'*Ex_bar; rle=eye(mu+1)*Ex_bar-rxy*inv(ryy)*rxy'; [v d]=eig(rle); d_temp=diag(d)'; [m,n ]= min(d_temp); b=norm_p*v(:,n)'; w=b*rxy*inv(ryy); conv(w,p); mmse=b*rle*b'; %by using easy way to get error energy c=conv(w,p) alpa=c(delta+1)./b(1) biased_error_energy=norm_p^2*(m-Ex_bar*(1-alpa)^2) unbiased_error_energy=norm_p^2*(m-Ex_bar*(1-alpa)^2)/alpa^2 SNRmfb=norm_p^2*Ex_bar/unbiased_error_energy; SNRmfb_in_dB=10*log10(SNRmfb) else %IIR %in case of IIR filter P matrix is infinite dimesion but there is a trick to get exact rxy, ryy norm_p=sqrt(1/(p(1)^2*(1-(p(2)/p(1))^2 ))); [temp size_p]=size(p); ptemp=[(-p(2)/p(1)).^(0:1:L-1)]/p(1); %-1 factor! P=toeplitz([ptemp(1) zeros(1,L-1)]',ptemp); Ptemp=toeplitz(ptemp',ptemp); ryy=Ex_bar*Ptemp*norm_p^2+sigma*eye(L,L); rxy=[zeros(mu+1,delta) eye(mu+1) zeros(mu+1,L-1-mu-delta)]*P'*Ex_bar; rle=eye(mu+1)*Ex_bar-rxy*inv(ryy)*rxy'; [v d]=eig(rle); d_temp=diag(d)'; [m,n ]= min(d_temp); b=norm_p*v(:,n)'; w=b*rxy*inv(ryy); c=conv(w,p); sum(conv(w,p).^2)-norm(b)^2*(w(1)/b(1))^2; mmse=b*rle*b'; %by using easy way to get error energy alpa=c(delta+1)./b(1) biased_error_energy=norm_p^2*(m-Ex_bar*(1-alpa)^2) unbiased_error_energy=norm_p^2*(m-Ex_bar*(1-alpa)^2)/alpa^2 SNRmfb=norm_p^2*Ex_bar/unbiased_error_energy; SNRmfb_in_dB=10*log10(SNRmfb); end