空間平滑MUSIC演算法（1）

1. 前言

在上一篇博客中有提到，當多個入射信號相干時，傳統MUSIC演算法的效果就會不理想，具體原因是多個入射信號相干時，有部分能量就會散發到噪聲子空間，使得MUSIC演算法不能對其進行有效估計，
針對這種情況，解相干方法被提出，本文主要講解通過降維處理解相干，之所以叫降維處理是因為這種方法將原陣列拆分成很多個子陣列，通過子陣列的協方差矩陣重構接收資料協方差矩陣，陣列的自由度DOF會隨之降低，即可分辨的相干信號數降低，
先看看傳統MUSIC演算法對相干信號進行DOA估計的效果，
MATLAB代碼

% Code For Music Algorithm
% The Signals Are All Coherent
% Author:癢羊羊
% Date：2020/10/29

clc; clear all; close all;
%% -------------------------initialization-------------------------
f = 500;                                        % frequency
c = 1500;                                       % speed sound
lambda = c/f;                                   % wavelength
d = lambda/2;                                   % array element spacing
M = 20;                                         % number of array elements
N = 100;                                        % number of snapshot
K = 6;                                          % number of sources
coef = [1; exp(1i*pi/6);... 
        exp(1i*pi/3); exp(1i*pi/2);... 
        exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals

%% generate signal
dd = (0:M-1)'*d;                                % distance between array elements and reference element
A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
X = A*(coef*S);                                 % received data without noise, M*N
X = awgn(X,10,'measured');                      % received data with SNR 10dB

%% calculate the covariance matrix of received data and do eigenvalue decomposition
Rxx = X*X'/N;                                   % covariance matrix
[U,V] = eig(Rxx);                               % eigenvalue decomposition
V = diag(V);                                    % vectorize eigenvalue matrix
[V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
U = U(:,idx);                                   % reset the eigenvector
P = sum(V);                                     % power of received data
P_cum = cumsum(V);                              % cumsum of V

%% define the noise space
J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
J = J(1);                                       % number of principal component
Un = U(:,J+1:end);

%% music for doa; seek the peek
theta = -90:0.1:90;                             % steer theta
doa_a = exp(-1i*2*pi*dd*sind(theta)/lambda);    % manifold array for seeking peak
music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
music = 10*log10(music/max(music));             % normalize the result and convert it to dB

%% plot
figure;
plot(theta, music, 'linewidth', 2);
title('Music Algorithm For Doa', 'fontsize', 16);
xlabel('Theta(°)', 'fontsize', 16);
ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
grid on;

可以看到，對于相干信號，傳統MUSIC演算法DOA估計效果很差，

2. 空間平滑演算法

降維處理解相干方法主要包括空間平滑類處理演算法，而空間平滑類處理演算法可以分成前向空間平滑演算法(FSS)、后向平滑演算法(BSS)、前后向平滑演算法(FBSS)，正如上面所說，這些演算法的估計效果很好，但是損失了陣列孔徑，導致可分辨的相干信號數降低，

2.1 線性陣列信號模型

關于這個線性陣列中符號的具體說明可以參考上一篇博文，子空間方法——MUSIC演算法

2.2 前向空間平滑演算法

前向空間平滑演算法將陣列分成多個互相重疊的子陣列，然后對子陣接收資料的協方差矩陣作平均，當子陣陣元數目大于等于相干信號數時，可以有效解相干，

如上圖所示，我們把M元陣列均勻分成L個子陣列，每個子陣列有N=M-L+1個陣元，以最左邊的子陣列為參考陣列，定義第J個子陣的接收資料為：
X J f = [ x J , x J + 1 , ? ? , x J + N ? 1 ] X_{J}^{f}=\left[x_{J}, x_{J+1}, \cdots, x_{J+N-1}\right] XJf?=[xJ?,xJ+1?,?,xJ+N?1?]那么第J個子陣接收資料的協方差矩陣(也稱作空間平滑矩陣)可以表示為
R N j j = A 1 D j ? 1 R s ( D j ? 1 ) H A 1 H + σ 2 I N j = 1 , 2 , ? ? , L R_{N}^{j j}=A_{1} D^{j-1} R_{s}\left(D^{j-1}\right)^{H} A_{1}^{H}+\sigma^{2} I_{N} \quad j=1,2, \cdots, L RNjj?=A1?Dj?1Rs?(Dj?1)HA1H?+σ2IN?j=1,2,?,L其中， D = diag ? [ e ? j ? 2 π d λ sin ? θ 1 , e ? j 2 π d λ sin ? θ 2 , ? e ? j 2 π d λ sin ? θ R ] D=\operatorname{diag}\left[\mathrm{e}^{-j\frac{-2 \pi d}{\lambda} \sin \theta_{1}}, \mathrm{e}^{-j\frac{2 \pi d}{\lambda} \sin \theta_{2}}, \cdots \mathrm{e}^{-j\frac{2 \pi d}{\lambda} \sin \theta_{R}}\right] D=diag[e?jλ?2πd?sinθ1?,e?jλ2πd?sinθ2?,?e?jλ2πd?sinθR?]，A1為第一個子陣即參考陣列的流型矩陣，
因此，前向空間平滑后的協方差矩陣可以由各個子陣的協方差矩陣求平均得到，
R ~ f = 1 L ∑ j = 1 L R N j j \tilde{R}^{f}=\frac{1}{L} \sum_{j=1}^{L} R_{N}^{j j} R~f=L1?j=1∑L?RNjj?利用前向空間平滑協方差矩陣和MUSIC演算法就可以分辨出多個相干信號的方位，可以證明，該方法最大可以檢測M/2個相干的信號，

MATLAB代碼

% Code For Music Algorithm
% The Signals Are All Coherent
% Author:癢羊羊
% Date：2020/10/29

clc; clear all; close all;
%% -------------------------initialization-------------------------
f = 500;                                        % frequency
c = 1500;                                       % speed sound
lambda = c/f;                                   % wavelength
d = lambda/2;                                   % array element spacing
M = 20;                                         % number of array elements
N = 100;                                        % number of snapshot
K = 6;                                          % number of sources
L = 10;                                         % number of subarray
L_N = M-L+1;                                    % number of array elements in each subarray
coef = [1; exp(1i*pi/6);... 
        exp(1i*pi/3); exp(1i*pi/2);... 
        exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals

%% generate signal
dd = (0:M-1)'*d;                                % distance between array elements and reference element
A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
X = A*(coef*S);                                 % received data without noise, M*N
X = awgn(X,10,'measured');                      % received data with SNR 10dB

%% reconstruct convariance matrix
%% calculate the covariance matrix of received data and do eigenvalue decomposition
Rxx = X*X'/N;                                   % origin covariance matrix
Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
for i = 1:L
    Rf = Rf+Rxx(i:i+L_N-1,i:i+L_N-1);
end
Rf = Rf/L;
[U,V] = eig(Rf);                                % eigenvalue decomposition
V = diag(V);                                    % vectorize eigenvalue matrix
[V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
U = U(:,idx);                                   % reset the eigenvector
P = sum(V);                                     % power of received data
P_cum = cumsum(V);                              % cumsum of V

%% define the noise space
J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
J = J(1);                                       % number of principal component
Un = U(:,J+1:end);

%% music for doa; seek the peek
dd1 = (0:L_N-1)'*d;
theta = -90:0.1:90;                             % steer theta
doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
music = 10*log10(music/max(music));             % normalize the result and convert it to dB

%% plot
figure;
plot(theta, music, 'linewidth', 2);
title('Music Algorithm For Doa', 'fontsize', 16);
xlabel('Theta(°)', 'fontsize', 16);
ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
grid on;

可以看到，在6個入射信號均相干的情況下，基于前向平滑的MUSIC演算法能較好地對其進行DOA估計，但仍存在估計精度的問題，比如真實入射角度為75°的信號的方位被估計成了74.2°，
后向空間平滑MUSIC演算法和前/后向空間平滑MUSIC演算法在下一篇博客中講解，

參考論文：《基于空間平滑MUSIC演算法的相干信號DOA估計，陳文鋒，吳桂清》
代碼Code均為原創，
歡迎轉載，表明出處，