《优化阵列信号处理上》鄢社锋-第二章代码

clc;clear;close all; c = 340; %声速 f = 1000; %频率 lambda = c/f; k = 2*pi/lambda; %波数 space = lambda/2; %阵元间距 M = 10; theta_d1 = 10*pi/180; theta_d2 = -10*pi/180; theta_angle = -90:0.1:90; theta = theta_angle*pi/180; p1 = exp(1i*k*space*[0:M-1]*sin(theta_d1)).'; p2 = exp(1i*k*space*[0:M-1]*sin(theta_d2)).'; Rx = 10^(30/10)*(p1*p1')+10^(0/10)*eye(M); %p1*p1'与p1'*p1不一样，这个地方卡了好久 %10° W_MVDR1 = Rx\p1 /(p1'*(Rx\p1)); B1 = zeros(size(theta)); for i = 1:length(theta) p = exp(1i*k*space*[0:M-1]*sin(theta(i))).'; B1(i) = W_MVDR1'*p; end B_dB1 = 20*log10(abs(B1)); %-10° W_MVDR2 = Rx\p2 /(p2'*(Rx\p2)); B2 = zeros(size(theta)); for i = 1:length(theta) p = exp(1i*k*space*[0:M-1]*sin(theta(i))).'; B2(i) = W_MVDR2'*p; end B_dB2 = 20*log10(abs(B2)); figure(1) plot(theta_angle,B_dB1,'LineWidth',1.5,'LineStyle','-.');hold on; plot(theta_angle,B_dB2,'LineWidth',1.5); grid on; xlabel('方位(^o)');ylabel('20lg(B)/dB');title('波束图');