clc;
clear;
r=0.5;                % space of array 
sensor=8;             % number of sensor
s_s=[r*(0:(sensor-1)/sensor);zeros(1,sensor)];
c=340;                % velocity of sound in air
t=10;                % 观测时间
fl=200;fh=300;        % 频带范围
f0=(fl+fh)/2;         % 中心频率
fs=1000;              % 采样频率
snap=fs*t;            % 采样点数
K=100;                 % 频域块拍
snap_s=snap/K;        % 频域点数
snaps=1000;           % fft点数

%----随机高斯信号通过带通滤波器--------
snr=10^(20/10);        % 信噪比
doa_azimuth=[20 0]*pi/180;                       % DOA
nu=length(doa_azimuth);                           % 信号个数
p=[cos(doa_azimuth);sin(doa_azimuth)];            % 方向向量
[b,a]=butter(5,[fl fh]/fs*2);                     % 带通滤波器
x=snr*(randn(nu,snap)+j*randn(nu,snap));          
for n=1:nu
    y(n,:)=filter(b,a,x(n,:));                    % 入射信号
end
df=(fh-fl)/snap_s;
fj=fl:4*df:fh;
k=length(fj)-1;
fl_s=fj(1:k);
fh_s=fj(2:k+1);
f0_s=(fh_s+fl_s)/2;

Y=0;
for u=1:k
    [b,a]=butter(5,[fl_s(u) fh_s(u)]/fs*2);
    for n=1:nu
        x_s(n,:)=filter(b,a,y(n,:));
    end
    A_s=exp(j*2*pi*f0_s(u)*s_s'*p/c);
    y_s=A_s*x_s;
    Y=Y+y_s;                       % 阵元接收到的信号   
end

Y=Y+randn(sensor,snap)+j*randn(sensor,snap);
f=(1:snaps)/snaps*fs;
f00=f0*1.1;                % 聚焦频率
f0_s=[f00,f0_s];

for k=1:length(fj)
    for l=1:length(f)
        if f(l)==f0_s(k);              % 所需频点对应位置
           q(k)=l;
        end
    end
end
YX=cell(1,K);
FYX=cell(1,K);
for n=1:K
YX{1,n}=Y(:,n:n+snap-K);             % 阵列采样输出分段
FYX{1,n}=fft(YX{1,n},snaps,2);       % 分段数据ｆｆｔ
Q{1,n}=FYX{1,n}(:,q);                % 提取频点数据

end
Q=cell2mat(Q);
z0=Q(:,1:k:k*(K-1)+1);
R0=z0*z0'/K;
[U0,S0,V0]=svd(R0);
Us0=U0(:,1:nu);                      % 聚焦频率对应的信号子空间

R=0;
for m=2:length(fj)
    z=Q(:,m:k:k*(K-1)+m);
    Rz=z*z'/K;
    [Uz,Sz,Vz]=svd(Rz);
    Usz=Uz(:,1:nu);
    T=Us0*Usz';                    % 聚焦矩阵
    Rt=T*Rz*T';
    R=R+Rt;
end
[U,S,V]=svd(R);
Us=U(:,1:nu);
Un=U(:,(nu+1):sensor);
theta=-90:0.5:90;
for n=1:length(theta)
    A=exp(j*2*pi*f00*s_s'*[cos(theta(n)*pi/180);sin(theta(n)*pi/180)]/c);
    P(n)=inv(A'*Un*Un'*A);
end
plot(theta, 20*log10(P))
这是我的程序！是用信号子空间法求的聚焦矩阵，请各位帮我找找问题所在吧 ！不胜感激！