SPA Matlab Code(转载)

最新推荐文章于 2022-03-13 11:01:38 发布
最新推荐文章于 2022-03-13 11:01:38 发布
阅读量3.8k 收藏 25
点赞数 1
分类专栏： MATLAB编程 文章标签： Matlab 代码转载 SPA 稀疏和参数算法

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

 一、均匀线阵主程序

% test of SPA in the ULA case

clear all
close all
clc


K = 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 case

clear all
close all
clc

K = 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 2013

if nargin < 3 || max(Omega) == length(Omega)  % ULA
    Omega = [];
else
    if length(Omega) ~= size(Y,1)
        error('Dimension of Omega not match!');
    end
end
if nargin < 2
    mode_noise = 2;
end

[M, N] = size(Y);

if N > 1
    Rhat = Y * Y' / N;
end

cvx_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);

end



function [u, sig] = SPA_nonsing(Rhat, mode_noise, Omega)
% SPA when Rhat is nonsingular

if isempty(Omega)         % ULA
    M = size(Rhat, 1);
    
    Rhatinv = inv(Rhat);
    TconjR = zeros(M,1);
    TconjR(1) = trace(Rhatinv);
    for j = 2:M
        TconjR(j) = 2*sum(diag(Rhatinv,j-1));
    end
    
    Rhathalf = sqrtm(Rhat);
    
    switch mode_noise
            
        case 1,  % equal noise variances
            cvx_solver sdpt3
            cvx_begin sdp            
              variable x(M,M) hermitian,
              variable u(M) complex,
              
              [x Rhathalf; Rhathalf toeplitz(u)] >= 0,
              
              minimize trace(x) + real(TconjR'*u);
            cvx_end
            sig = 0;
            
        case 2,  % different noise variances
            cvx_solver sdpt3
            cvx_begin sdp
              variable 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_end
            
        otherwise,
            error('error!');
    end
    
    return
end


%% SLA

M = 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:M
    TconjR(j) = 2*sum(diag(Rhatinv_expand,j-1));
end

S = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);

Rhathalf = sqrtm(Rhat);

switch mode_noise
        
    case 1,  % equal noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable 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_end
        
        sig = 0;
        
    case 2,  % different noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable 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_end
        
    otherwise,
        error('error!');
end

end




function [u, sig] = SPA_SMV(Y, mode_noise, Omega)
% SPA in the single-snapshot case

if isempty(Omega)         % ULA
    M = length(Y);
    
    Y2n = norm(Y);
    
    switch mode_noise
            
        case 1,  % equal noise variances
            cvx_solver sdpt3
            cvx_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_end
            
            sig = 0;
            
        case 2,  % different noise variances
            cvx_solver sdpt3
            cvx_begin sdp
              variable 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_end
            
        otherwise,
            error('error!');
    end
    
    return
end


%% SLA

M = max(Omega);
Mbar = length(Omega);

Y2n = norm(Y);

S = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);

switch mode_noise
        
    case 1,  % equal noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable x
          variable u(M) complex,
          
          toeplitz(u) >= 0,
          [x Y2n*Y'; Y2n*Y S*toeplitz(u)*S'] >= 0,
            
          minimize x + Mbar*real(u(1));
        cvx_end
        
        sig = 0;
        
    case 2,  % different noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable x
          variable 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_end
        
    otherwise,
        error('error!');
end

end







function [u, sig] = SPA_singu(Rhat, mode_noise, Omega)
% SPA when Rhat is singular

if isempty(Omega)         % ULA
    M = size(Rhat,1);
    
    switch mode_noise
            
        case 1,  % equal noise variances
            cvx_solver sdpt3
            cvx_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_end
            
            sig = 0;
            
        case 2,  % different noise variances
            cvx_solver sdpt3
            cvx_begin sdp
              variable 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_end
            
        otherwise,
            error('error!');
    end
    
    return
end


%% SLA

M = max(Omega);
Mbar = length(Omega);

S = zeros(Mbar, M);
S(:, Omega) = eye(Mbar);

switch mode_noise
        
    case 1,  % equal noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable 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_end
        
        sig = 0;
        
    case 2,  % different noise variances
        cvx_solver sdpt3
        cvx_begin sdp
          variable 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_end
        
    otherwise,
        error('error!');
end

end




function [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 amp

prec = 1e-4; % can be tuned appropriately

eigval = sort(real(eig(toeplitz(u))));
if eigval(1) > 0
    sigma_in_u = eigval(1);
    u(1) = u(1) - sigma_in_u;
    eigval = eigval - sigma_in_u;
else
    sigma_in_u = 0;
end

K = 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

 

ArcheL
关注
专栏目录
特征波长优选SPA特征matlab代码.zip
05-12
1.版本：matlab2014/2019a/2021a 2.附赠案例数据可直接运行matlab程序。 3.代码特点：参数化编程、参数可方便更改、代码编程思路清晰、注释明细。 4.适用对象：计算机，电子信息工程、数学等专业的大学生课程设计、期末大作业和毕业设计。
SPA matlab程序
11-01
SPA程序，可以运行，有数据，很直观，可与GUI的说明配合使用！！！
端元提取——逐次投影算法SPA的MATLAB代码
04-29
自己学习时编写的混合像元分解中的端元提取部分的逐次投影算法SPA的MATLAB代码，浅显易懂。可以运行得到结果，但没有结果用于验证，亦未编写精度评定，不保证结果的正确性。
【ELM分类】基于SPA特征极限学习机实现建筑材料光谱分类附matlab代码
matlab_dingdang的博客
03-13 1252
​1 简介 近年来,随着多媒体技术、数字化等技术的快速发展,机器学习领域得到了迅速的发展。极限学习机(ELM)作为一个有竞争力的机器学习算法,以其简洁的理论和易于实施的特点吸引了越来越多学者的关注。在分类和回归等问题中,极限学习机得到了广泛的应用,并表现出了良好的效果。 2 部分代码 %% 基于SPA+ELMclcclear allclose all%% 读取数据% 红砖data1=load('红砖.txt');data1=data1';% % 石灰% data2=load('石灰.tx
spa.zip_livinglfp_matlab的SPA_plsr_spa matlab_spa的matlab实现
07-15
使用matlab，spa算法提取光谱的特征波长，结合了plsr
matlab连续投影算法SPA使用示例
热门推荐
PeanutbutterBoh的博客
07-04 2万+
连续投影算法（Successive Projections Algorithm，SPA）是一种使矢量空间共线性最小化的前向变量选择算法， 它的优势在于提取全波段的几个特征波长，能够消除原始光谱矩阵中冗余的信息，可用于光谱特征波长的筛选。
【BP分类】基于SPA特征BP神经网络实现光谱分类附matlab代码
matlab_dingdang的博客
03-13 1462
1 简介 BP神经网络用于光谱分类的具体流程如图1所示,其步骤具体如下: 1)采集样本数据集并对数据做预处理,通常BP网络的输入数据都需要做归一化处理。 2)将已知数据样本分为训练集和检验集两部分,训练样本集用于对神经网络进行训练,检验样本集用于对训练好的网络进行检验。 3)选择神经网络的结构和规模,主要是选择网络的层数和每层的单元数。 4)选择一定的学习算法用训练样本集对网络进行训练,使网络的实际输出与理想输出满足一定的误差范围。 5)用检验样本集对网络进行检...
spa matlab
最新发布
07-29
- *2* *3* [SPA Matlab Code(转载)](https://blog.csdn.net/qq_25634581/article/details/85683399)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_...
极化码极化现象matlab代码.rar_polar_polar code_信道极化 现象_极化码 极化_极化码代码
07-14
2. **构造极化矩阵**：极化码的核心是极化矩阵（BCH矩阵），它是由一系列的比特串并行-串行转换（BSC）和串行-并行转换（SPA）构成的。在MATLAB中，可以使用循环或递归方式构建这个矩阵。 3. **信息位到编码位的...
LDPC码编译码matlab代码
06-27
**LDPC码（Low-Density Parity-Check Code）**是一种高效的纠错编码技术，由Robert G. Gallager在1962年提出，但真正受到广泛关注是在20世纪90年代，由于其接近香农极限的性能而备受瞩目。在通信、数据存储等领域...
MATLAB uve建模 cars spa波段挑选程序
07-28
MATLAB uve建模 cars spa波段挑选程序 用于光谱或数据的分类和回归
高斯白噪声matlab代码-SPA_for_LDPC:此存储库是关于LDPC（又名低密度奇偶校验）代码的和积算法在二进制对称信道，二进制擦除信
05-26
高斯白噪声matlab代码SPA_for_LDPC 这个存储库是关于LDPC（又名低密度奇偶校验）代码的和积算法（在二进制对称信道，二进制擦除信道和AWGN（加性高斯白噪声）下）的实现（使用C和Matlab）的） 渠道。 感谢您在中提供这些（几乎）常规LDPC矩阵文件。 感谢Takuji Nishimura和devoloping The，也感谢Shawn Cokus提供了。
spa
02-14
自述文件 该自述文件通常会记录启动和运行应用程序所需的所有步骤。 您可能要讲的内容： Ruby版本 系统依赖 配置 数据库创建 数据库初始化 如何运行测试套件 服务（作业队列，缓存服务器，搜索引擎等） 部署说明 ...
SPA-连续投影算法-教程版
04-11
对老版本代码进行了调试，找了一个光谱的数据进行降维测试，效果不错，有需要的小伙伴下载吧，调试不易，谢谢支持。
matlab的compan,MATLAB经常使用操做
weixin_33402252的博客
03-19 1740
1.点乘，点除，点乘方javascript点乘(对应元素相乘),必须同维或者其中一个是标量，a.*bjava点除，a.\b表示矩阵b的每一个元素除以a中对应元素或者除以常数a，a./b表示常数a除以矩阵b中每一个元素或者矩阵a除以矩阵b对应元素或者常数b算法点乘方a.^b,矩阵a中每一个元素按b中对应元素乘方或者b是常数数组2.矩阵中元素的操做函数矩阵a中第r行，a(r,:),spa第r列，a(:...
opencv读取matlab标定,相机标定与矫正opencv+MATLAB
weixin_42651748的博客
03-22 708
本文目的在于记录如何使用MATLAB作摄像机标定,并经过opencv进行校订后的显示。ios首先关于校订的基本知识经过OpenCV官网的介绍便可简单了解：ide对于摄像机咱们所关心的主要参数为摄像机内参，以及几个畸变系数。上面的链接中后半部分也给了如何标定，然而OpenCV自带的标定程序稍显繁琐。于是在本文中我主推使用MATLAB的工具箱。下面让咱们开始标定过程。spa标定板方法二：逼格满满(MA...
【SVM分类】基于SPA特征支持向量机SVM实现光谱分类附matlab代码
matlab_dingdang的博客
03-13 2394
​1 简介 目前基于紫外-可见光光谱法的建筑物废弃材料判别模型大多使用主成分分析(PCA)进行特征提取,然而,对于光谱相似度较高的有建筑物废弃材料,仅从数据驱动的角度提取其方差最大的方向作为特征进行识别效果往往不佳。针对有机污染物光谱数据多重共线性以及谱峰重叠干扰的问题,开展了基于连续投影算法(SPA)和支持向量机(SVM)的紫外-可见光光谱建筑物废弃材料判别方法研究。实验结果表明, SPA作为一种提取光谱数据原始特征波段的方法,可以有效的对建筑物废弃材料的紫外-可见光光谱进行特征提取,提升不同物质之间的
花一天时间搜索,居然找到一个德国人写的论文里面有用Matlab读取*.spa文件的程序片段
stereohomology
05-22 3152
虽然是德语的, 但是, 1.程序大部分都是英文; 2.剩余的小量的描述, 用google翻译搞定就可以了 http://att.newsmth.net/att.php?p.455.56499.365.png 有待于用Matlab实际检验之...
spafdr matlab,MATLAB在微阵列数据分析的FDR控制中的实现
weixin_34853926的博客
03-16 375
0、引言在多重假设检验问题中,首先是将多个单重假设检验作为一个整体(称为一个检验族),然后从总体上对这个检验族中的每个假设同时进行检验。狭义的多重比较是对多组的总体参数或多组的分布比较后的组间两两比较,比如方差分析后多组均值的比较等。而多重检验是多重比较的推广。多重比较一般指多变量情形,即对同一问题通过对多变量的反复检验来回答,如微阵列数据分析中各个因子的组间比较的假设检验等,多元回归中各自变量的...
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值