%% 阵列信号处理第一次作业
clc;close all;clear all;
j = sqrt(-1);
fs = 1000; %系统采样频率
T = 0.4;
t = 0:1/fs:0.4; %信号持续时间
M = 16; %阵元个数
[nn,N] = size(t);
f = 62.5; %入射信号频率
c = 340; %声速
lamda = c / f; %信号波长
d = lamda / 2; %阵元间距
ss = sin(2*pi*f*t);
sc = cos(2*pi*f*t);
s = ss+j*sc; %信号表示
w = 1/M * ones(M,1); %均匀加权的归一化取值
%% 问题(a)
Angle = 0:10:100; %信号入射角度范围
[nn,n] = size(Angle); %得到不同角度个数
ya = zeros(n,N); %频率-波束响应空间
i = 1;
for ang = Angle
fi = 2*pi*d*cos(ang*pi/180)/lamda; %波数Kz与d的乘积
v = exp(j*fi*(0:M-1)'); %线列阵的频率-波束响应
w = 1/M * ones(M,1); %均匀加权的归一化取值
yy =w' * (v * s); %均匀加权线列阵的输出
ya(i,:) = yy - mean(yy);
i = i + 1;
end
figure
i = 1:11;
plot(t,ya(i,:));
legend('0度','10度','20度','30度','40度','50度','60度','70度','80度','90度','100度');
title(['不同入射角度下波束形成器的输出序列波形']);
xlabel('t/s');ylabel('幅度');
ylim([-1.2 1.2]);
figure
for i = 1:4
subplot(4,1,i)
plot(t,ya(i,:));title([num2str((i-1)*10),'度方向时的输出序列']);xlabel('t');ylabel('y');
end
figure
for i = 5:8
subplot(4,1,i-4)
plot(t,ya(i,:));title([num2str((i-1)*10),'度方向时的输出序列']);xlabel('t');ylabel('y');
end
figure
for i = 9:11
subplot(3,1,i-8)
plot(t,ya(i,:));title([num2str((i-1)*10),'度方向时的输出序列']);xlabel('t');ylabel('y');
end
%% 问题(b)
a = 0;
for SNR = 10:-10:-10 %信噪比的三个取值
Angle2 = 0:180; %信号入射角度范围
[m,n] = size(Angle2); %得到不同角度个数
yb = zeros(1,n);
i = 1;
for ang2 = Angle2
fi = 2*pi*d*cos(ang2*pi/180)/lamda; %波数Kz与d的乘积
v = exp(j*fi*(0:M-1)'); %线列阵的频率-波束响应
x = zeros(M,N);
for k = 1:M
xx = v(k) * s; %线列阵的输出信号
x(k,:) = xx;
end
yb(i) = w' * x * (w' * x )'; %均匀加权线列阵的输出信号功率
i = i + 1;
end
yb = 10*log10(yb/max(yb)); %对输出分贝归一化处理
plot(Angle2,yb); %绘输出分贝与角度的关系
title('无噪声时的分贝输出');xlabel('角度');ylabel('分贝');axis([0 180 -50 0]);
end
%% 问题(c)
a = 0;
for SNR = 10:-10:-10 %信噪比的三个取值
Angle2 = 0:180; %信号入射角度范围
[m,n] = size(Angle2); %得到不同角度个数
yb = zeros(1,n);
k = 1;
for ang = Angle2
tao = abs(cos(ang*pi/180) )/2/f;
x = zeros(M,N);
for i = 1:M
temp = (i-1)*tao;
xx = sin(2*pi*f*(t-temp)) .* (stepfun(t,temp) - stepfun(t,temp+T));
x(i,:) = awgn(xx,SNR);
end
yb(k) = w' * x * (w' * x )'; %均匀加权线列阵的输出信号功率
k = k + 1;
end
yb=10*log10(yb/max(yb)); %对输出分贝归一化处理
a = a + 1;
subplot(3,1,a);
plot(Angle2,yb); %绘输出分贝与角度的关系
title(['输出序列平方求和后的分贝输出SNR=',num2str(10-(a-1)*10),'dB']);xlabel('角度');ylabel('分贝');
end
%% 问题(d)
R = 20; %旁瓣比主波束低26dB
x0 = cosh(1/(M-1)*acosh(R)); %求出最大值对应的x轴值
fip = 2 * acos( (1/x0) * cos((2*(1:M-1)-1)*pi/(M-1)/2) );%得到零点
FI = [0 fip]; %阵列流形矩阵
V = exp(j * (0:M-1)' * FI);
e1 = zeros(M,1);
e1(1) = 1;
w = inv(V') * e1; %权值矢量
a = 0;
figure
for SNR = 10:-10:-10; %信噪比的三个取值
Angle3 = 0:180; %信号入射角度范围
[m,n] = size(Angle3); %得到不同角度个数
yc = zeros(1,n);
k = 1;
for ang = Angle3
x = zeros(M,N);
tao = abs(cos(ang*pi/180) )/2/f;
for i = 1:M
temp = (i-1)*tao;
xx = sin(2*pi*f*(t-temp)) .* (stepfun(t,temp) - stepfun(t,temp+T));
x(i,:) = awgn(xx,SNR);
end
yc(k) = w' * x * (w' * x )'; %均匀加权线列阵的输出信号功率
k =k + 1;
end
yc=10*log10(yc/max(yc)); %对输出分贝归一化处理
a = a + 1;
subplot(3,1,a);
plot(Angle3,yc); %绘输出分贝与角度的关系
title(['DC加权输出序列平方求和后的分贝数输出',num2str(SNR),'dB']);xlabel('角度');ylabel('分贝');
end
%% 问题(e)
theta = 0:180; %信号扫描角度范围
for k = 1:length(theta) %取角度大小
for n=0:M-1
a(k,n+1) = exp(-1j*(n-(M-1)/2)*2*pi*d/lamda*cos(theta(k)*pi/180));
end
end
w = 1/M; %均匀加权的归一化取值
s = exp(1j*2*pi*d/lamda*cos(60*pi/180)*[0:M-1]')*exp(1j*2*pi*f*t);
y = w*a*s;
result = y*y';
figure;
out = diag(result);
plot(theta,10*log10(out/max(out)));
title('无噪声时入射信号60°方向的均匀加权波束形成结果');
xlabel('\theta/°');
ylabel('幅度');
axis([0 180 -50 0]);
%% 问题(f)
w = 1/M; %均匀加权的归一化取值
signal_direction = 60;
for SNR = -10:10:0; %信噪比的两个取值
s = exp(1j*2*pi*d/lamda*cos(60*pi/180)*[0:M-1]')*exp(1j*2*pi*f*t);
s = awgn(s,SNR);
y = w*a*s;
result = y*y';
figure
out = diag(result);
plot(theta,10*log10(out/max(out)));
title(['SNR=',num2str(SNR),'时入射信号60°方向的均匀加权波束形成结果']);
xlabel('\theta/°');
ylabel('幅度');
end
评论0