Matlab实现空间信息和极化信息估计的MUSIC算法

clear; close all; clc; M=6; signum=2; N=1000; theta1=80*pi/180; theta111=25; %0.4363 [0,2*pi] theta2=40*pi/180; theta222=40; % 0.6981 fai1=50*pi/180; fai111=50; % 0.8727 [0,pi] fai2=60*pi/180; fai222=60; %1.0472 gama1=70*pi/180; gama111=35; %0.6109 [0,pi/2] gama2=20*pi/180; gama222=30; %0.5236 eta1=60*pi/180; jita111=15; % 0.2618[-pi,pi] eta2=50*pi/180; jita222=20; % 0.3491 f1=0.30; f2=0.50; snr=10; MM=1; %----s--- s1=44*(exp(j*2*pi*(f1*0.1*[0:N-1]+rand))); % fc/fs=0.7*0.1 s2=44*(exp(j*2*pi*(f2*0.1*[0:N-1]+rand))); % fc/fs=1.0*0.1 s=[s1(1:N);s2(1:N)]; % a11=v_a(theta111,fai111,gama111,jita111); % a22=v_a(theta222,fai222,gama222,jita222); a1=[-sin(theta1),cos(theta1)*cos(fai1);... cos(theta1),sin(theta1)*cos(fai1);... 0,-sin(fai1);... cos(theta1)*cos(fai1),sin(theta1);... sin(theta1)*cos(fai1),-cos(theta1);... -sin(fai1),0]*[cos(gama1);sin(gama1)*exp(j*eta1)]; a2=[-sin(theta2),cos(theta2)*cos(fai2);... cos(theta2),sin(theta2)*cos(fai2);... 0,-sin(fai2);... cos(theta2)*cos(fai2),sin(theta2);... sin(theta2)*cos(fai2),-cos(theta2);... -sin(fai2),0]*[cos(gama2);sin(gama2)*exp(j*eta2)]; A =[a1,a2]; noise=randn(6,N)+j*randn(6,N); z=awgn(A*s,snr); for hang=1:6, for lie=1:6, zmat = z(1,:).* conj(z(1,:)).* z(hang,:); R4 = (zmat * z(lie,:)') /N; r12= (z(1,:) * z(1,:)') / N; r34= (z(hang,:)*z(lie,:)')/N; r13= (z(1,:) * z(hang,:).') / N; r24= (conj(z(1,:))*z(lie,:)')/N; r14= (z(1,:) * z(lie,:)') / N; r23= (conj(z(1,:))*z(hang,:).')/N; C(hang,lie) = R4 - r12*r34-r13*r24-r14*r23; end end [u s v]=svd(C); En=v(3:6,:)'; En=u(:,3:6); step=2; theta=(0:step:360)*(pi/180); phi=(0:step:180)*(pi/180); Pmusic=zeros(length(theta),length(phi)); for i=1:length(theta) for j=1:length(phi) D=[-sin(theta(i)),cos(theta(i))*cos(phi(j));... cos(theta(i)),sin(theta(i))*cos(phi(j));... 0,-sin(phi(j));... cos(theta(i))*cos(phi(j)),sin(theta(i));... sin(theta(i))*cos(phi(j)),-cos(theta(i));... -sin(phi(j)),0]; Pmusic(i,j)=1/(det(D'*(En)*En'*D)); end end Pmusic=abs(Pmusic); Pmusic=Pmusic/max(max(Pmusic)); Pmusic=10*log10(Pmusic); figure(1); mesh(phi*(180/pi),theta*(180/pi),Pmusic); %colormap([0 0 0]); xlabel('\phi'); ylabel('\theta'); zlabel('Nomalized Power Spectrum of Music/dB'); figure(2); contour(phi*(180/pi),theta*(180/pi),Pmusic),grid on; %colormap([0 0 0]); xlabel('\phi'); ylabel('\theta'); gama=(0:step:90)*(pi/180); eta=(-180:step:180)*(pi/180); % theta1=theta1+0.01; % theta2=theta2+0.02; Pmusic=zeros(length(gama),length(eta)); D1=[-sin(theta1),cos(theta1)*cos(fai1);... cos(theta1),sin(theta1)*cos(fai1);... 0,-sin(fai1);... cos(theta1)*cos(fai1),sin(theta1);... sin(theta1)*cos(fai1),-cos(theta1);... -sin(fai1),0]; D2=[-sin(theta2),cos(theta2)*cos(fai2);... cos(theta2),sin(theta2)*cos(fai2);... 0,-sin(fai2);... cos(theta2)*cos(fai2),sin(theta2);... sin(theta2)*cos(fai2),-cos(theta2);... -sin(fai2),0]; for i=1:length(gama) for j=1:length(eta) gg=[cos(gama(i));sin(gama(i))*exp(1i*eta(j))]; Z1=D1'*(En)*En'*D1; Pmusic1(i,j)=1/(det(gg'*D1'*(En)*En'*D1*gg)); Z2=D2'*(En)*En'*D2; Pmusic2(i,j)=1/(det(gg'*D2'*(En)*En'*D2*gg)); end end Pmusic1=abs(Pmusic1); Pmusic1=Pmusic1/max(max(Pmusic1)); Pmusic1=10*log10(Pmusic1); Pmusic2=abs(Pmusic2); Pmusic2=Pmusic2/max(max(Pmusic2)); Pmusic2=10*log10(Pmusic2); figure(3); mesh(eta*(180/pi),gama*(180/pi),Pmusic1); %colormap([0 0 0]); xlabel('\eta'); ylabel('\gamma'); zlabel('Nomalized Power Spectrum of Music/dB'); figure(4); contour(eta*(180/pi),gama*(180/pi),Pmusic1); grid on; %colormap([0 0 0]); xlabel('\eta'); ylabel('\gamma'); figure(5); mesh(eta*(180/pi),gama*(180/pi),Pmusic2); %colormap([0 0 0]); xlabel('\eta'); ylabel('\gamma'); zlabel('Nomalized Power Spectrum of Music/dB'); figure(6); contour(eta*(180/pi),gama*(180/pi),Pmusic2); grid on; %colormap([0 0 0]); xlabel('\eta'); ylabel('\gamma'); % gama_11 = % 69.9948 % gama_22 = % 20.0056 % eta_11 = % 40.0514 % eta_22 = % 60.0334