matlab_数字移动通信中子带滤波器组处理

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % demo6_7_QP2bf.m % K. Bell 12/2/99 % updated by K. Bell 11/20/00 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Orthogonal QP constraints %************************ % Array %************************ N = 10; % Elements in array d = 0.5; % sensor spacing half wavelength wrt wc D = [-(N-1)/2:1:(N-1)/2].'; BWNN = 2/(N*d); u=[-1:0.001:1]; us = [0 [0.001:0.001:0.01] [0.015:0.005:0.1]]; ns = length(us); nu=length(u); vv = exp(j*2*pi*d*D*u); uI = 0.3; vI = exp(j*2*pi*d*D*uI); sigma_i = 10^(20/10); sigma_n = 1; sigma_s = 10^(-20/10); DL = 10^(20/10); DL=0; % Dolph Chebychev SL = -25; wt = poly(exp(j*2*acos(cos((2*[1:1:N-1]-1)*pi/(2*(N-1)))/cosh(acosh(10^(-SL/20))/(N-1))))).'; wdq = real(wt/sum(wt)); u_d = 0.1; %wdq = u_d*sinc(D*u_d); wdqbar = wdq/real(wdq'*wdq); Pq_orth = eye(N)-wdq*inv(wdq'*wdq)*wdq'; r = 2*u_d*sinc(-[0:1:N-1]*u_d); c = 2*u_d*sinc([0:1:N-1]*u_d); Q = toeplitz(c,r); R = Pq_orth*Q*Pq_orth; [v,lam] = eig(R); [lam,ind] = sort(diag(abs(lam))); % sort eigenvalues in ascending order v = v(:,ind); % arrange eigenvectors in same order C = [wdqbar v(:,N) v(:,N-1)]; % QP +first 2 three eigenvectors g = [1 0 0].'; %************************ % Source %************************ AS = exp(j*2*pi*d*D*us); SINRin = sigma_s/(sigma_i+sigma_n); T = zeros(1,ns); Tc = zeros(1,ns); Tn = zeros(1,ns); I = zeros(1,ns); Ic = zeros(1,ns); In = zeros(1,ns); A = zeros(1,ns); An = zeros(1,ns); Ac = zeros(1,ns); O = zeros(1,ns); On = zeros(1,ns); Oc = zeros(1,ns); figure(1) whitebg('k') figure(2) whitebg('k') for n=1:ns Rs = sigma_s*AS(:,n)*AS(:,n)'; wc = AS(:,1)/N; Bc = 10*log10(abs(wc'*vv).^2); Snc = sigma_i*vI*vI'; Sn = Snc+sigma_n*eye(N); Sx = Rs+Sn; Sxinv = inv(Sx+DL*eye(N)); Sninv = inv(Sn+DL*eye(N)); w = Sninv*C*inv(C'*Sninv*C)*g; wn = Sxinv*C*inv(C'*Sxinv*C)*g; T(n) = 10*log10(real(w'*w)); Tc(n) = 10*log10(real(wc'*wc)); Tn(n) = 10*log10(real(wn'*wn)); O(n) = 10*log10(real(w'*Rs*w)); Oc(n) = 10*log10(real(wc'*Rs*wc)); On(n) = 10*log10(real(wn'*Rs*wn)); I(n) = 10*log10(real(w'*Snc*w)); Ic(n) = 10*log10(real(wc'*Snc*wc)); In(n) = 10*log10(real(wn'*Snc*wn)); A(n) = 10*log10(real((w'*Rs*w)/(w'*Sn*w))/SINRin); Ac(n) = 10*log10(real((wc'*Rs*wc)/(wc'*Sn*wc))/SINRin); An(n) = 10*log10(real((wn'*Rs*wn)/(wn'*Sn*wn))/SINRin); B = 10*log10(abs(w'*vv).^2); Bn = 10*log10(abs(wn'*vv).^2); figure(1) h=plot(u,B); hold on hc=plot(u,Bc,'g'); set(hc,'Linewidth',1.0) hp=plot(u,Bn,'m'); set(hp,'Linewidth',1.0) hh=plot(us(:,n)*[1 1],[-60 30],'c'); set(hh,'Linewidth',1.0) hh=plot(uI*[1 1],[-60 30],'r'); set(hh,'Linewidth',1.0) xlabel('u') ylabel('Beampattern (dB)') title(['QP2 Const., Signal Power ' num2str(10*log10(sigma_s)) ' dB, DL = ' num2str(10*log10(DL)) ' dB']) grid on hold off axis([-1 1 -50 20]) legend([hc h hp],'Conventional','LCMV','LCMP') figure(2) subplot(4,1,1) h1=plot(us(1:n),A([1:n]),'y'); set(h1,'Linewidth',1.0) hold on h2=plot(us(1:n),Ac([1:n]),'g'); set(h2,'Linewidth',1.0) h3=plot(us(1:n),An([1:n]),'m'); set(h3,'Linewidth',1.0) hold off axis([0 0.1 -30 30]) grid on xlabel('u_a') ylabel('Array gain (dB)') legend([h2 h1 h3],'Conventional','LCMV','LCMP',3) title(['QP2 Const., Signal Power ' num2str(10*log10(sigma_s)) ' dB, DL = ' num2str(10*log10(DL)) ' dB']) subplot(4,1,2) h=plot(us(1:n),O([1:n]),'y'); set(h,'Linewidth',1.0) hold on h=plot(us(1:n),Oc([1:n]),'g'); set(h,'Linewidth',1.0) h=plot(us(1:n),On([1:n]),'m'); set(h,'Linewidth',1.0) hold off SDB = 10*log10(sigma_s); v=[0 0.1 SDB-50 SDB+10]; axis(v) grid on xlabel('u_a') ylabel('Signal Pwr. (dB)') subplot(4,1,3) h=plot(us(1:n),T([1:n]),'y'); set(h,'Linewidth',1.0) hold on h=plot(us(1:n),Tc([1:n]),'g'); set(h,'Linewidth',1.0) h=plot(us(1:n),Tn([1:n]),'m'); set(h,'Linewidth',1.0) hold off axis([0 0.1 -20 20]) grid on xlabel('u_a') ylabel('Sensit. (Noise Pwr.)(dB)') subplot(4,1,4) h=plot(us(1:n),I([1:n]),'y'); set(h,'Linewidth',1.0) hold on; h=plot(us(1:n),Ic([1:n]),'g'); set(h,'Linewidth',1.0) h=plot(us(1:n),In([1:n]),'m'); set(h,'Linewidth',1.0) hold off axis([0 0.1 -80 20]) grid on xlabel('u_a') ylabel('Interference Pwr. (dB)') drawnow %pause end