Music算法仿真

Music：利用信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，估计信号频率。这里的谱峰也称为“伪谱”，或MUSIC谱。

以下列信号为例  

 

• 产负高斯白噪声序列和带噪声的信号样本

N=1000；
noise=(rand(1,N)+j*rand(1,N))/sqrt(2);
u(n)=exp(1i*0.5*pi*(0:N-1)+1i*2*pi*rand)+exp(-1i*0.3*pi*(0:N-1)+1i*2*pi*rand)+noise;

 

• 计算自相关矩阵

M=8;%自相关矩阵的阶数
for k=1:N-M
   xs(:,k)=un(k+M-1:-1:k).';
end
R=xs*xs'/(N-M);

• 计算伪谱（考虑到要画点，就要记录每次P的离散值，因此w的范围用n来表示，下面用NF个采样点表示）

NF=2048;
for n=1:NF
    aw=exp(-1i*2*pi*((n-1)/(NF-1)-0.5)*(0:M-1)');
    P(n)=1/(aw'*G*G'*aw);
    f(n)=((n-1)/(NF-1)-0.5);
end
maxx=max(P);
figure,plot(-f,10*log10((P+eps)/maxx));
xlabel('频率f');ylabel('归一化功率谱/dB')  

注意这里绘图时f要取负号，因为当最后个点n=NF时，(n-1)/(NF-1)-0.5 相当于f=-0.5，但是我们画的坐标轴f=-0.5在x轴最左边，因此要取个符号。

   

 

 

完整代码

%% 产生 1000 个样本数的随机过程
clear all,close all
N=1000;%样本数
M=8;K=2;
noise=(rand(1,N)+1i*rand(1,N))/sqrt(2);
n=1:1:N;
u(n)=exp(1i*0.5*pi*n+1i*2*pi*rand(1))+exp(-1i*0.3*pi*n+1i*2*pi*rand(1))+noise;
%% 计算自相关矩阵（计算时可用 FFT 计算）
u2=[u,zeros(1,N)]; %将 u(n)进行扩展，补 N 个零
U2=fft(u2);
P=abs(U2).^2/N; %计算功率谱
r2=ifft(P);
r=r2(1:M);%取前 M 个值
rc=r2(2*N-M+2:2*N);%取出 r(-M+1):r(-1)
rc=[rc,r(1)];%r(-M+1):r(0)
rc=fliplr(rc);%首尾颠倒
Rxx=toeplitz(r,rc)';
%% 对 Rxx 进行特征分解
[V,D]=eig(Rxx);%特征值为什么没有虚部，还是为 complex
D=abs(diag(D));%因为我看到 D 是复数，虽然没有虚部，
[Dnew,I] = sort(D); %对其排序，小到大排列
Vnew = V(:, I) ;
G=Vnew(:,1:M-K); %M-K 个噪声子空间的特征向量， G 即为该空间
%% Music 算法估计
NF=2048;
for n=1:NF
aw=exp(-1i*2*pi*((n-1)/(NF-1)-0.5)*(0:M-1)'); %n=NF 时，相当于此时 f=-0.5
P(n)=1/(aw'*G*G'*aw);
f(n)=((n-1)/(NF-1)-0.5);
end
maxx=max(P);
figure,plot(-f,10*log10((P+eps)/maxx));
xlabel('频率 f');ylabel('归一化功率谱/dB')
%% 求解赋初值
syms z;
pz = z.^([0:M-1]');
pz1 = (z^(-1)).^([0:M-1]);
fz = z.^(M-1)*pz1*G*G'*pz;
% fz=aN*G*G'*aP';如果是这个，在下一步会报错说不是多项式，因为含有 z 的负数次幂
factor=sym2poly(fz);
zx=roots(factor);
rx=zx.';
[as,ad]=(sort(abs((abs(rx)-1))));
frequency=angle(rx(ad([1,3])))./(2*pi)