%RSS算法(聚焦变量为阵列流型矩阵) 可修改阵元数和阵元间距
clc
clear all
close all
%% 常数定义
deg2rad = pi/180; % deg -> rad
rad2deg = 180/pi; % rad -> deg
N_array = 30; % 阵元数
N_snap = 2000; % 快拍数
fc=100; %中心频率
fl = 80; % 入射LFM信号最低频率
fu = 120; % 入射LFM信号最高频率
fs = 4*fu; % 采样频率(说明:这里设置的fu和fs恰好能保证信号在相邻阵元间的延迟为一个快拍点)
T = 1/fs; % 采样时间间隔
t = T:T:N_snap*T;% 采样时间点
Ts=N_snap*T; %信号持续时间
K=(fu-fl)/Ts;
snr=50; %信噪比
c = 3e8; % 声速
lamda = c/fu; % 波长
d = lamda/4; % 阵元间距
ang_deg = [80 120]; % 信源方位角(角度)
ang_rad = ang_deg * deg2rad; % 信源方位角(弧度)
P=length(ang_deg); %信源数
%kk=20; %时域分段数
%s(1,:)=exp(ima*2*pi*(fc*t+0.5*K*t.^2)); % lfm信号
%s(2,:)=exp(ima*2*pi*(fc*t-0.5*K*t.^2));
s(1,:)=exp(j*2*pi*(fl*t+0.5*K*t.^2)); % lfm信号1
s(2,:)=exp(j*2*pi*(fu*t+0.5*K*t.^2)); %lfm信号2
%% 阵列接收数据
%duan=N_snap/kk; %时域每一段包含duan(100)个点
data = zeros(N_array, length(t));%阵列接收的数据
signal = zeros(N_array, length(t));
for n=1:N_array
signal(n,:) = 1*exp(1j*2*pi*(fl+1/2*K*t).*t).*exp(1j*2*pi./(c./(fl+1/2*K*t))*d.*(n-1)*sin(ang_rad(1)))+1*exp(1j*2*pi*(fu+1/2*K*t).*t).*exp(1j*2*pi./(c./(fu+1/2*K*t))*d.*(n-1)*sin(ang_rad(2)));%假设原点处信号的相位为0;%假设原点处信号的相位为0
end
data=signal;
%SS=awgn(data,00,'measured'); %加白噪声
n=exp(j*2*pi*randn(N_array,N_snap))/snr;%噪声
SS=data+n;
SS_FFT = zeros(N_array,N_snap*N_snap);
for seg=1:N_snap %数据段
for m = 1:N_array
SS_FFT(m,N_snap*(seg-1)+1:N_snap*seg)=fft(SS(m,1*(seg-1)+1),N_snap); %时域不分段 共2000个时间点 对每个点做2000次fft变换 产生2000*2000个频率数据 存入SS_FFT中
end
end
%% 选取聚焦矩阵
fy=repmat(fc,[1 P]);
for k=1:N_array
for sig=1:P
A0(k,sig)=exp(j*2*pi*fy(sig)*(k-1)*d*sin(ang_rad(sig))/c); % 确定聚焦频率的阵列流型矩阵 30*2(阵元数*信源数)
end
end;
%ii=0;
fXX=zeros(N_array,20000);
for fn=300:10:N_snap/4 %频率范围72Hz~120Hz
f=fs/(N_snap)*fn; %FFT后的频率 频率分辨率0.24Hz
f2=fs/(N_snap)*(fn+5); %子带中心频率
for i=1:N_snap
fXX=SS_FFT(:,fn+N_snap*(i-1):(fn+10)+N_snap*(i-1));
end
% fXX=SS_FFT(:,(fn+5):N_snap:(N_snap-1)*N_snap+(fn+5)); %对每个频点建立类似时域的阵列接收数据 30*20(阵元数*时域分段数)
% RfXX=fXX*fXX'/size(fXX,2);%频域的协方差阵 30*30(阵元数*阵元数)
% for jj=1:8
%f=fl+(jj-1)*5;
%fn=f*N_snap/fs;
% fXX=SS_FFT(:, fn:N_snap:(N_snap-1)*N_snap+fn);
RfXX=fXX*fXX'/size(fXX,2);%频域的协方差阵
Y=zeros(N_array,N_array);
for sig=1:P
%A(:,sig)=exp(j*2*pi*fff(jj)*(0:M-1)'*d*cos(angle(sig)*pi/180)/c);
A(:,sig)=exp(j*2*pi*f2*[0:N_array-1]'*d*cos(ang_rad(sig))/c); %聚焦变量是阵列流型 30*2(阵元数*信源数)
end
AA=A*A0'; % 30*30(阵元数*阵元数)
[U S V]=svd(AA);
%求确定频点的聚焦矩阵
T=V*U';
Y=Y+T*RfXX*T';
end
%% MUSIC算法
[EigenVectors,EigenValues]=eig(Y);
Lemda=diag(EigenValues); %计算矩阵特征值
[SortedLemda,Index]=sort(Lemda);
Index=flipud(Index); %将特征值降序排列
NoiseSubspace(1:N_array,1:N_array-P)=EigenVectors(1:N_array,Index(P+1:N_array));%把噪声空间分成两部分,分别为加入信号空间的一部分和噪声
%%在0到90之间进行搜索
for i=1:1:180
a=exp(j*2*pi*fc*d*[0:N_array-1]'*cos(i*pi/180)/c);
MUSIC_Spec(i)=a'*a/(a'*NoiseSubspace*(NoiseSubspace)'*a);
end
M_MUSIC_Spec=max(MUSIC_Spec);%把最大值选出
abslog_MUSIC_Spec=20*log10(abs(MUSIC_Spec));%/max(abs(MUSIC_Spec)
%谱峰搜索
ang=[];
index=find(abslog_MUSIC_Spec<60);%谱峰搜索的门限
abslog_MUSIC_Spec(index)=0;
figure
MUSIC=abs(MUSIC_Spec)/max(abs(MUSIC_Spec));
Delta=[1:180];
plot(Delta,20*log10(abs(MUSIC_Spec)),'r'); %NonNormative Spatial Spectrum
hold on;
%stem(ang,abslog_MUSIC_Spec(ang))
xlabel('角度(度)'),ylabel('空间方位谱(DB)');
grid on;
tic;
两个lfm宽带信号,频率范围80~120,时域不分段,共2000个时间点,对每个时间点做2000次快速傅里叶变换,聚焦矩阵聚焦在中心频率100HZ上,频率分辨率0.24Hz,十个点作为一段,子带2.4Hz,使用RSS方法聚焦。
找不到哪里有问题,谢谢