Root-Music方法进行频率估计

Root-Music算法是在Music算法的基础上提出的。

Music算法:利用信号子空间和噪声子空间的正交性,构造空间谱函数,通过谱峰搜索,估计信号频率。

具体就是先对信号用矩阵形式表示出来,接着求解相关矩阵,并对其进行特征分解,得出有K个非零特征值,利用剩下的M-K个零特征值得出信号子空间与噪声子空间的正交性。而实际中,自相关矩阵R是估计出来的,因此特征向量本身也带有一定误差,因此最终的公式中有求和(平滑作用)。

Root-Music算法:直接构造函数,并将噪声子空间的向量写成矩阵G,最终将信号频率估计问题转化成了一元高次方程的求根问题。在实际求解时,需要在2(M-1)个根中,找出其中位置最接近单位元的K个根,其相位就是信号频率的估计值。


  • FFT计算自相关函数:补零,FFT,计算功率谱,IFFT得II,然后截取前面m个值N个即为我们的信号自相关函数r(0),r(1)...r(N)
  • 计算相关矩阵:构造M阶相关矩阵来说,对于复值信号来说II的最后M-1即为r(-7)~r(-1),用fliplr()左右镜像为r(-1)~r(-7),II的前M个即为r(0)~r(7),然后就是用toeplitz来计算。。(网上看到好几个用xcorr直接计算相关系数,然后用toeplitz)
%% 计算自相关矩阵(计算时可用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)';
  • 对相关矩阵进行特征值分解:因为上面求出的相关矩阵是一个复数类型,matlab调用eig求解后特征值为复数类型(相关矩阵的特征值为非负实数,看了下里面的数都是虚部为0i,然而matlab却还是认为是复数),难道每次都要对特征值取abs吗?      QAQ
  • 得到M-K个特征子空间,即对特征值进行排列,取出最小的M-K个特征值和特征向量。(eigs可解决?)
    if ~issorted(diag(D)) [V,D] = eig(A); [D,I] = sort(diag(D)); V = V(:, I) 
  • 接下来我在想如何把方程写出来,因为这是一个矩阵G,因此变量也得表示成一个向量。。
    syms z
    az = z.^([0:M-1]');
  • 最后对整个表达式进行系数提取(注意的是如果直接提取会报错sym2poly输入非多项式,因为这是含有z的负数次幂,因此可等价转化乘以最低次),并求解,求出来是复数解,找出其中最接近单位圆的K个根(最接近单位圆可能有很多根,最后选的时候稍微隔开几个),这些根的相位angle函数就是信号频率的估计。
    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])))./pi


  • 3
    点赞
  • 44
    收藏
    觉得还不错? 一键收藏
  • 6
    评论
以下是一个简单的二维root-music频率估计的MATLAB代码: ```matlab clc; clear; % 生成数据 N = 500; % 采样点数 M = 5; % 信号数 theta = [30 40 50 60 70]; % 信号角度 fs = 1000; % 采样频率 T = 1/fs; % 采样时间间隔 L = 10; % 阵元数 d = 0.5; % 阵元间距 lambda = 1; % 波长 Omega = 2*pi*d/lambda; % 阵元间距的相位差 % 构造阵列 array = zeros(L, 2); for ii = 1:L array(ii,:) = [(ii-1)*d 0]; end % 生成信号 t = (0:N-1)*T; s = zeros(N, M); for ii = 1:M s(:,ii) = exp(1j*2*pi*theta(ii)*cos(deg2rad(array(:,1)))'*sin(deg2rad(array(:,2))))'; end % 加入噪声 SNR = 10; % 信噪比 noise = (randn(N,L) + 1j*randn(N,L))/sqrt(2); noise_power = norm(s(:))/sqrt(N*M*10^(SNR/10)); % 噪声功率 x = s + noise*noise_power; % 二维root-music估计 Rxx = x'*x/N; [EV,D] = eig(Rxx); [y,i] = sort(diag(D),'descend'); V = EV(:,i(M+1:end)); theta_range = -90:0.1:90; % 角度范围 Pmusic = zeros(length(theta_range),length(theta_range)); for ii = 1:length(theta_range) for jj = 1:length(theta_range) a = exp(1j*2*pi*Omega*d*cos(deg2rad(theta_range(ii)))*(0:L-1)'); b = exp(1j*2*pi*Omega*d*sin(deg2rad(theta_range(jj)))*(0:L-1)'); Pmusic(ii,jj) = 1/(a'*V*V'*b); end end % 绘制估计结果 figure(1); subplot(1,2,1); plot(array(:,1),array(:,2),'o','MarkerSize',10,'LineWidth',2); axis([-2 2 -2 2]); axis square; xlabel('x'); ylabel('y'); title('阵列'); subplot(1,2,2); contour(theta_range,theta_range,10*log10(abs(Pmusic).^2),[-80:5:0]); colormap jet; colorbar; xlabel('\theta_x (degree)'); ylabel('\theta_y (degree)'); title('二维root-music估计结果'); ``` 代码的主要思路是生成一个二维阵列,通过在阵列上采集信号,然后加入噪声,最后使用二维root-music算法进行频率估计。代码中的主要步骤包括: 1. 生成阵列:使用给定的阵元数和阵元间距,生成一个二维阵列。 2. 生成信号:根据给定的信号角度、采样频率和阵列,生成一个采样点数为N、信号数为M的信号矩阵。 3. 加入噪声:根据给定的信噪比,生成与信号矩阵大小相同的噪声矩阵,并将其与信号矩阵相加。 4. 二维root-music估计:根据估计的角度范围,计算每个角度的似然函数值,并使用等高线图展示结果。 代码中使用了MATLAB自带的contour函数绘制等高线图。你可以通过修改代码中的参数来测试不同情况下的频率估计效果。
评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值