目录
1. 根据正弦信号参数估计确定单目标测向系统中的CRLB
1.1 实数正弦信号参数估计问题
根据 [1] 例题3.14中关于正弦信号参数估计的问题，如果采集的数据满足：
其中信号幅度A，频率fs ,起始相位φ 都为未知参数，N 为采样点数。
根据该例题及（3.31）中的结论，可以得到：
根据信噪比的定义，为信号功率和噪声功率的比值，在归一化的负载上，对于实数信号的功率可以表示为
,而噪声信号功率为
。因此实数信号时，存在：
因此，
也可以改写为：
1.2 复正弦信号参数估计问题
如果采集的信号不是实数信号，而是复数信号，即：
根据[1] 例题15.13中关于复正弦信号参数估计问题，可以得到：
在归一化的负载上，对于复数信号的功率可以表示为
,而噪声信号功率为
。因此复数信号时，存在：
所以，对复正弦信号进行参数估计，估计频率的CRLB满足：
由于后续讨论都采用复数信号，因此直接以复正弦信号的CRLB用于讨论，并记
，那么得到参数都未知情况下频率的CRLB为：
需要注意的是，上述信噪比不是对数格式（dB）的信噪比。
1.3?利用上述结论计算传统单目标测向系统的CRLB
图1?单目标线性天线阵列测向模型?
转变到测向问题时，模型如上图所示。假设当前仅有唯一测向目标，天线1获取的信号为：
那么天线2~M获取的信号可以近似表示为：
……
其中：
可以看成由于接收天线空间分别有关的相位延时，θ 为入射信号方向，l 为天线间距, f0 为测向信号频率。
根据《统计信号处理基础-估计与检测理论》例题3.15中的处理方式，假设了每根天线采集一个快拍数据，那么总共能够获得M个数据，分别对于的相位为：
0，
，……,?
?
因此例题3.15中的方向估计问题可以等效为15.13中的频率估计问题，此时3.15中的
与15.13中的
存在如下关系：
因此，由参数变化后的CRLB性质：
上式为测向模型中，由于位置引起的相位延时的CRLB下限。
如果需要求信号入射方向的CRLB，此时3.15中的
与15.13中的
存在如下关系：
也就是：
因此：
也就是：
考虑到线性阵列的总长度为：
因此，上式可以变为：
上式为测向模型中，估计唯一入射信号方向的CRLB下限。
2. 根据阵列测向模型确定单目标测向系统中的CRLB
多目标测向时，线性阵列的模型可以描述为：
其中
那么根据已证明的结论[2]-[3]，均匀线阵的CRLB：
其中，N是每个天线的采样点数，
如果是单目标时,且每个天线采样1个点，即N=1
此时：
又由于：
?那么可以得到：
而又因为信号功率可以表示为：
此时信噪比也可以表示为：
因此，单目标，且每个天线单快拍采样情况下，
的CRLB下限可以表示为：
这个结果与用正弦信号方式推导的测向问题一致。因此对于单目标测系统，两种方法都可以使用。最后，对于单目标测向系统，如果考虑采样点数，那么信号延时和入射方向的CRLB下限可以分别表示为：
3.?Matlab仿真示例及结果
示例代码如下，每种信噪比，仿真了10000次求方差
图2?单目标，单快拍，4天线线阵情况下MUSIC测向在不同信噪比下的方差与CRLB对比，注意，横坐标为SNR为dB情况下的值
图3?将上述计算结果去lg对数，可以获得进行线性的递增效果?
?通过对比，可以发现Music测向的结果，在目前仿真的信噪比范围上，能近似接近CRLB下限，也表明了改方法的性能。
clc;clear; close all;
%%%%%%%% MUSIC CRLB for Uniform Linear Array%%%%%%%%
derad = pi/180; %角度->弧度
M = 4; % 阵元个数
K = 1; % 信源数目
N = 1; % 快拍数
delay=0.1*pi; %待估计的延时，可以由测向信号产生
snr=10:2:30; %遍历信噪比范围
trail=10000; %尝试次数
d=0:1:(M-1); %线阵间隔
esti_music=zeros(1,trail);
Len_SNR=length(snr);
VAREsti_snr=zeros(Len_SNR,1);
VARCalc_snr=zeros(Len_SNR,1);
for index_snr=1:Len_SNR
snr_current=snr(index_snr)
A=exp(-j*d.'*delay); %方向矢量，复数形式
for index_trail=1:trail
S=randn(K,N); %信源信号，入射信号，不相干即可，也可以用正弦替代
X=A*S; %构造接收信号
X1=awgn(X,snr_current,'measured'); %引入高斯白噪声，此时的SNR为对数形式
Rxx=X1*X1'/N; %标准MUSIC算法，计算协方差矩阵
[EV,D]=eig(Rxx); %特征值分解
EVA=diag(D)'; %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA); %将特征值排序 从小到大
EV=fliplr(EV(:,I)); % 对应特征矢量排序
iang = -pi:0.01*pi:pi; %延时遍历范围
for index=1:length(iang)
angle_input=iang(index);
phim=angle_input;
a=exp(-j*d*phim).';
En=EV(:,K+1:M); % 取矩阵的第M+1到N列组成噪声子空间
Pmusic(index)=1/(a'*En*En'*a); %行程MUSIC谱
end
Pmusic=abs(Pmusic);
[y x]=max(Pmusic); %单目标，求最大值即可
esti_music(index_trail)=iang(x(1)); %寻找最大值对应延时
end
VAREsti_snr(index_snr)=(sum((esti_music-delay).^2)/(trail)); %计算方差
SNR_Linear=10^(snr_current/10); %将对数SNR转换成线性格式
VARCalc_snr(index_snr)=6/(N*M*(M-1)*(M+1)*SNR_Linear);
end
figure(1) %线性CRLB作图，近似为1/x函数
plot(snr,VAREsti_snr,'b')
hold on
plot(snr,VARCalc_snr,'r')
hold off
figure(2) %对数CRLB作图，近似斜率为负的线性函数
plot(snr,10*log(VAREsti_snr),'b')
hold on
plot(snr,10*log(VARCalc_snr),'r')
hold off
?参考文献
1.?Steven?M.Kay, 统计信号处理基础-估计与检测理论，2014
2.?刘永良，空间谱估计理论及算法， 2004
3.?P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and Cramer-Rao bound," in?IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 5, pp. 720-741, May 1989
