SPA Matlab Code(转载)

注释：代码转载自 https://sites.google.com/site/zaiyang0248/publication

一、均匀线阵主程序

% test of SPA in the ULA caseclear all
close all
clcK = 3;          % source number
M = 30;         % array length
N = 200;        % snapshot number
sigma = 1;      % noise variance% true frequency and power
f = [.1; .11; .5];
p = [3; 3; 1];A = exp(1i*2*pi*kron((0:M-1)',f'));
S = sqrt(diag(p))*exp(1i*2*pi*rand(K,N));
% S(K,:) = S(1,:); % sources 1 & K coherent
E = diag(sqrt(sigma/2)) * (randn(M,N)+1i*randn(M,N));% observed snapshots
Y = A*S + E;% parameter estimation using SPA
tic
mode = 1;  % noise variance mode, 1: equal; 2: different
[freq_est, pow_est] = SPA(Y, mode);
toc% MSE
MSE = norm(f - sort(freq_est(1:K)))^2 / K;% plot result
figure(100), semilogy(f, p, 'ko'), hold on;
semilogy(freq_est, pow_est, 'rx', 'MarkerSize', 5);
xlabel('Frequency');
ylabel('Power');
legend('Truth', 'Estimate');
title(strcat('Frequency MSE = ', num2str(MSE)));

二、稀疏线阵主程序

% test of SPA in the SLA caseclear all
close all
clcK = 6;
Omega = [1; 2; 5; 7];
M = max(Omega);
L = length(Omega);
T = 200;% true frequency and power
f = [.1; .18; .4; .55; .7; .85];
p = [2; 2; 2; 1; 1; 1];% observed snapshots
X = exp(1i*2*pi*kron(Omega-1,f'));
S = diag(sqrt(p))*exp(1i*2*pi*rand(K,T));sigma = .1;
E = sqrt(sigma/2) * (randn(L,T)+1i*randn(L,T));
Y = X*S + E;% parameter estimation using SPA
tic
mode = 1;  % noise variance mode, 1: equal; 2: different
[freq_est, pow_est] = SPA(Y, mode, Omega);
toc% MSE
MSE = norm(f - sort(freq_est(1:K)))^2 / K;% plot result
figure(100), plot(f, p, 'ko'), hold on;
plot(freq_est, pow_est, 'rx', 'MarkerSize', 5);
xlabel('Frequency');
ylabel('Power');
% axis([0,1,0,3]);
legend('Truth', 'Estimate');
title(strcat('Frequency MSE = ', num2str(MSE)));

三、SPA函数

function [freq, pow] = SPA(Y, mode_noise, Omega)% [freq, pow] = SPA(Y, mode_noise, Omega)
%
% SPA implements the sparse and parametric approach (SPA) for linear array
% signal processing.
%
% Input
% Y: observation
% mode_noise: 1 if equal noise variances;
%             2 if different noise variances,
%             default value: 2
% Omega: index set of SLA, valid only in the SLA case, sorted ascendingly
%
% Output:
% freq: frequency estimate
% pow: power estimate
% sigma: noise variance estimate
%
% Reference: Z. Yang, L. Xie, and C. Zhang, "A discretization-free sparse
%   and parametric approach for linear array signal processing",
%   IEEE Trans. Signal Processing, 2014
%
% Written by Zai Yang, April 2013if nargin < 3 || max(Omega) == length(Omega)  % ULAOmega = [];
elseif length(Omega) ~= size(Y,1)error('Dimension of Omega not match!');end
end
if nargin < 2mode_noise = 2;
end[M, N] = size(Y);if N > 1Rhat = Y * Y' / N;
endcvx_quiet true
cvx_precision default% solve the SDP (or estimate R)
if N >= M && cond(Rhat) < 1e8       % nonsingular[u, sig] = SPA_nonsing(Rhat, mode_noise, Omega);elseif N == 1  % single snapshot[u, sig] = SPA_SMV(Y, mode_noise, Omega);else           % multisnapshots and singular[u, sig] = SPA_singu(Rhat, mode_noise, Omega);end% postprocessing and parameter estimation
[freq, pow] = PostProc(u);endfunction [u, sig] = SPA_nonsing(Rhat, mode_noise, Omega)
% SPA when Rhat is nonsingularif isempty(Omega)         % ULAM = size(Rhat, 1);Rhatinv = inv(Rhat);TconjR = zeros(M,1);TconjR(1) = trace(Rhatinv);for j = 2:MTconjR(j) = 2*sum(diag(Rhatinv,j-1));endRhathalf = sqrtm(Rhat);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdp            variable x(M,M) hermitian,variable u(M) complex,[x Rhathalf; Rhathalf toeplitz(u)] >= 0,minimize trace(x) + real(TconjR'*u);cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable x(M,M) hermitian,variable u(M) complex,variable sig(M), sig >= 0,toeplitz(u) >= 0,[x Rhathalf; Rhathalf toeplitz(u)+diag(sig)] >= 0,minimize trace(x) + real(TconjR'*u) + real(diag(Rhatinv)')*sig;cvx_endotherwise,error('error!');endreturn
end%% SLAM = max(Omega);
Mbar = length(Omega);Rhatinv = inv(Rhat);
Rhatinv_expand = zeros(M);
Rhatinv_expand(Omega, Omega) = Rhatinv;
TconjR = zeros(M,1);
TconjR(1) = trace(Rhatinv_expand);
for j = 2:MTconjR(j) = 2*sum(diag(Rhatinv_expand,j-1));
endS = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);Rhathalf = sqrtm(Rhat);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdpvariable x(Mbar,Mbar) hermitian,variable u(M) complex,toeplitz(u) >= 0,[x Rhathalf; Rhathalf S*toeplitz(u)*S'] >= 0,minimize trace(x) + real(TconjR'*u);cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable x(Mbar,Mbar) hermitian,variable u(M) complex,variable sig(Mbar), sig >= 0,toeplitz(u) >= 0,[x Rhathalf; Rhathalf S*toeplitz(u)*S'+diag(sig)] >= 0,minimize trace(x) + real(TconjR'*u) + real(diag(Rhatinv)')*sig;cvx_endotherwise,error('error!');
endendfunction [u, sig] = SPA_SMV(Y, mode_noise, Omega)
% SPA in the single-snapshot caseif isempty(Omega)         % ULAM = length(Y);Y2n = norm(Y);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdp            variable x,variable u(M) complex,toeplitz(u) >= 0,[x Y2n*Y'; Y2n*Y toeplitz(u)] >= 0,minimize x + M*real(u(1));cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable x,variable u(M) complex,variable sig(M), sig >= 0,toeplitz(u) >= 0,[x Y2n*Y'; Y2n*Y toeplitz(u)+diag(sig)] >= 0,minimize x + M*real(u(1)) + sum(sig);cvx_endotherwise,error('error!');endreturn
end%% SLAM = max(Omega);
Mbar = length(Omega);Y2n = norm(Y);S = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdpvariable xvariable u(M) complex,toeplitz(u) >= 0,[x Y2n*Y'; Y2n*Y S*toeplitz(u)*S'] >= 0,minimize x + Mbar*real(u(1));cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable xvariable u(M) complex,variable sig(Mbar), sig >= 0,toeplitz(u) >= 0,[x Y2n*Y'; Y2n*Y S*toeplitz(u)*S'+diag(sig)] >= 0,minimize x + Mbar*real(u(1)) + sum(sig);cvx_endotherwise,error('error!');
endendfunction [u, sig] = SPA_singu(Rhat, mode_noise, Omega)
% SPA when Rhat is singularif isempty(Omega)         % ULAM = size(Rhat,1);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdp            variable x(M,M) hermitian,variable u(M) complex,toeplitz(u) >= 0,[x Rhat; Rhat toeplitz(u)] >= 0,minimize trace(x) + M*real(u(1));cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable x(M,M) hermitian,variable u(M) complex,variable sig(M), sig >= 0,toeplitz(u) >= 0,[x Rhat; Rhat toeplitz(u)+diag(sig)] >= 0,minimize trace(x) + M*real(u(1)) + sum(sig);cvx_endotherwise,error('error!');endreturn
end%% SLAM = max(Omega);
Mbar = length(Omega);S = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);switch mode_noisecase 1,  % equal noise variancescvx_solver sdpt3cvx_begin sdpvariable x(Mbar,Mbar) hermitian,variable u(M) complex,toeplitz(u) >= 0,[x Rhat; Rhat S*toeplitz(u)*S'] >= 0,minimize trace(x) + Mbar*real(u(1));cvx_endsig = 0;case 2,  % different noise variancescvx_solver sdpt3cvx_begin sdpvariable x(Mbar,Mbar) hermitian,variable u(M) complex,variable sig(Mbar), sig >= 0,toeplitz(u) >= 0,[x Rhat; Rhat S*toeplitz(u)*S'+diag(sig)] >= 0,minimize trace(x) + Mbar*real(u(1)) + sum(sig);cvx_endotherwise,error('error!');
endendfunction [freq, amp] = PostProc(u)% [freq, amp] = PostProc(u)
% PostProc implements the postprocessing procedure given the solution u.
% The result is sorted in the descending order of ampprec = 1e-4; % can be tuned appropriatelyeigval = sort(real(eig(toeplitz(u))));
if eigval(1) > 0sigma_in_u = eigval(1);u(1) = u(1) - sigma_in_u;eigval = eigval - sigma_in_u;
elsesigma_in_u = 0;
endK = sum(eigval > prec*eigval(end));M = length(u);% solve h
h = toeplitz([conj(u(2)); u(1:M-1)], conj(u(2:K+1))) \ (-u(1:M));% solve the zeros of H(z) in the form of exp(-1i * 2 * pi * theta)
r = roots([1; h]);
% r = r ./ abs(r);% solve the amper vector
mat = flipud(vander([r; zeros(M-K,1)]).');
amp = mat(:, 1:K) \ u;
[amp, idx] = sort(real(amp), 'descend');
Khat = sum(amp>0);
amp = amp(1:Khat);% determine the frequency
freq = - phase(r(idx(1:Khat))) / (2*pi);
freq = mod(freq, 1);end

SPA Matlab Code(转载)相关推荐

心电信号的PQRST模拟matlab代码(转载+自己调研汇总)
目前PQRST网上现成的代码有两份[1][2] 我们采用[2],原因是[1]中采用了lowpass这个函数, [3]中提到:"注意,只有2018年之后的matlab才有lowpass, ba ...
事件触发控制 Event-Trigger Control Matlab Code
提示:仅供参考事件触发控制Matlab Code 前言一.线性系统二.有限差分法三. Matlab ODE45 总结前言自动控制的大部分工作和研究考虑了周期或时间触发控制系统,在这些系统中 ...
期货策略matlab,code 一个利用MATLAB编写的螺纹钢期货高频交易套利策略联合开发网 - pudn.com...
code 所属分类:金融证券系统开发工具:matlab 文件大小:506KB 下载次数:398 上传日期:2013-10-09 14:14:53 上传者:huangxiao 说明: 一个利用M ...
matlab 高频数据,高频股价数据算波动率的matlab code (Yacine Ait-Sahalia)
用 R 算高频数据的日波动率可以用 package " highfrequency ",这里给出的是matlab code. 说明:高频数据的波动率code 这个文件夹对应的是A ...
ADC 动态参数分析matlab code的几个问题（span,spanh取值问题，幅度比例因子的添加等）
前言目前内外网自用ADC动态参数计算的matlab程序都是同一套模板的改版,存在共同的问题,迟迟没有得到解决,也困扰了我几天,不过还好,最终搞明白了几点. 一.ADC动态参数分析的源程序从何而来? ...
MATLAB Code of Artificial Potencial Field Method for Robot Path Planning 人工势场法局部极小问题
APF_Code download 模拟退火法处理局部极小问题源代码 MATLAB Code of Artificial Potencial Field Method for Robot Path P ...
如何让matlab提速,[转载]matlab提速技巧（自matlab帮助文件）
1.首先要学会用profiler.1.1. 打开profiler. To open the Profiler, select View -> Profiler from the MATLAB d ...
matlab速度梯度,[转载]关于FLUENT中Y+的一些讨论
一.关于 fluent计算时壁面函数法和网格的关系,还有一个小问题 1:各位用 fluent的同仁和高手们,我想要比较好的使用 fluent软件,最重要的就是要学好理论,在这里我想请教各位一个问题, ...
matlab profile,[转载]matlab profile的使用【原创】
用profile可解决的问题: 1.避免不必要的计算 2.改动代码避免耗时的函数 3.储存一定的结果,避免重复计算 profile用作调试工具: 1.查找出没有实际运行的代码 2.You can al ...

SPA Matlab Code(转载)

SPA Matlab Code(转载)相关推荐

最新文章

热门文章