窄带/宽带MUSIC方位估计

窄带MUSIC

宽带MUSIC

代码如下：

theta_s = 50;%	声源方位角
phi_s = 90;%	声源俯仰角
B = zeros( 91, 361 );%	方位谱
M = 16;%	阵元数
c = 1500;%	水下声速
f = 100;%	声波频率
fs = 10000;%	采样率
T = 1;
t = 0:1/fs:T;
lamda = c / f;
r = 2*lamda;%	圆形阵列半径
pos = r*[cos(2*pi*[1:M]/M); sin(2*pi*[1:M]/M); zeros(1, M) ];
k = 2*pi/lamda*[sind(phi_s)*cosd(theta_s); sind(phi_s)*sind(theta_s); cosd(phi_s)];
p = exp( 1i*k'*pos );
band = 0;
K = band/T;
SourceNum = 1;
% signal = p.'*exp( 1i*2*pi*f*t );%   cw
signal = p.'*exp( 1i*(2*pi*f*t+pi*K*t.^2) );%   lfm
N = 1;%	将原始信号分为N段

theta = linspace( 0 , 360 , 361 );% 方位角
phi = linspace( 0 , 90 , 91 );% 仰角
[m, n] = size( signal );
%--------------------------------------------------------------------------
%	窄带
% R = signal * signal' / n;
% % 特征值分解
% [ Q , D ] = eig( R );
% [ ~ , I ] = sort( diag( D ) , 1 , 'descend' );
% Q = Q( :,I );
% Qn = Q( : ,SourceNum + 1:end );                                  % 信号子空间
% for ii = 1:361
%     for jj = 1:91
%         k = 2 * pi / lamda * [ sind(phi(jj))*cosd(theta(ii));sind(phi(jj))*sind(theta(ii));cosd(phi(jj))];
%         p = exp( 1i*k'*pos );p = p.';
%         B( jj , ii ) = (1/(p'*(Qn*Qn')*p));
%     end
% end
% B = abs( B );
%--------------------------------------------------------------------------
%	宽带MUSIC
signalLength = ceil( n / N );
if mod( n, N ) ~= 0
    signal = [ signal, zeros( m, N*signalLength-n ) ];
end
signalFFT = zeros( m, 1024 );
for num1 = 1:N
    signalTemp = signal( :, (( num1-1 )*signalLength+1: num1*signalLength) );
    signalTempFFT = fft( signalTemp, 1024, 2 );
    signalFFT = signalFFT + signalTempFFT;
end
signalFFT = signalFFT / N;
fre = (0:1024/2-1)/1024*fs;
for num2 = 1:20:1024/2
    R = signalFFT( :, num2 )*signalFFT( :, num2 )';
    [ Q, D ] = eig( R );
    [ ~, I ] = sort( diag( D ), 1, 'descend' );
    Q = Q( :,I );
    Qn = Q( : ,SourceNum + 1:end );% 信号子空间
    FTemp = zeros( 91, 361 );
    lamda = c/fre(num2);
    for ii = 1:361
        for jj = 1:91
            k = 2 * pi / lamda * [ sind(phi(jj))*cosd(theta(ii)); sind(phi(jj))*sind(theta(ii)); cosd(phi(jj)) ];
            p = exp( 1i*k'*pos );p = p.';
            FTemp( jj, ii ) = (1/(p'*(Qn*Qn')*p));
        end
    end
    B = B + FTemp;
end
B = abs( B ) / signalLength;
%--------------------------------------------------------------------------

figure(1)
mesh(B)
xlabel( 'Azimuth angle/°' ); ylabel( 'Pitch angle/°' ); zlabel( 'Amp' );
set( gca ,'FontName','Times New Roman' );