例题3.6

代码

Main.m

clear all;clc;

N=256;M=4;K=3;

SNR1=30;SNR2=30;SNR3=27;%设定信噪比；

A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3);%求得信号的幅度；

f1=0.1;f2=0.25;f3=0.27;%设定频率值

s1=zeros(1,N);s2=zeros(1,N);s3=zeros(1,N);%初始化sk

s1=getsk(A1,f1,N);%求sk

s2=getsk(A2,f2,N);

s3=getsk(A3,f3,N);

v=randn(1,N);% 构建高斯白噪声；

x=s1+s2+s3+v;%构建了时域上的信号

N3=2048;%设定画图时描点的数目。

R=getR(x,N,M);

[G,D]=getG(R,M,K);

d=1/(N3-1);%求画图用的横坐标的间隔。

h=zeros(1,N3);

for i=1:N3

    h(i)=-0.5+(i-1)*d;

end

y=zeros(1,N3);

for j=1:N3

    w=h(j)*2*pi;

    aw=getaw(w,M);MinValue=min(y(j));

    y(j)=10*log10(abs(1/((aw')*G*(G')*aw)));

end

plot(h,y);

getAk.m 

function AK=getAk(SNR) %定义信号幅度

AK=((10^(SNR/10))*2)^0.5;

Getaw.m    %求aw

function aw=getaw(w,M)

aw=zeros(M,1);

for j=1:M

    aw(j)=exp(-w*(j-1)*i);

end

getR.m  %求自相关矩阵R

function R=getR(x,N,M)

L=N-M+1;

tempx=zeros(M,1);

R=zeros(M,M);

for n=M:N

    for j=1:M

        tempx(j)=x(n-(j-1));

    end

   R=R+tempx*tempx';

end

R=R/L;

getG.m  %求G矩阵

function [G,D]=getG(R,M,K)

[V,D]=eig(R);

G=zeros(M,M-K);

z=zeros(M,1);

for j=1:M

    z(j)=D(j,j);%将特征值放入了z里面

end

[z,y]=sort(z);%对z进行了排序目的是，找到最小的M-K个特征值多对应的特征向量。

for i=1:M

    for j=1:(M-K)

        G(i,j)=V((y(j)-1)*M+i);%第j小的特征值，对应的特征向量，是原来V中的第y(j)个列向量。

    end

end

getsk.m  %定义信号sk

function sk=getsk(Ak,fk,N)

qk=zeros(1,N);

% temp=rand(1,N);

% qk=temp*6.28*0.01;%设定随机相位的值

for n=1:N

    sk(n)=Ak*exp(6.28*fk*n*i+qk(n));

end

仿真图：

M=4

改变M：

M=8时

M=16时

由图可得M取值越大，得到的图像上的峰值对应的频率越精确。M=16时，MUSIC谱的峰值几乎精确地出现在f1=0.1,f2=0.25,f3=0.27,具有很高的分辨率。下载本文