等效源法（ESM，也称波叠加法）Matlab程序分享
等效源法简介：
所谓等效源法(Equivalent Source Method)，属于近场声全息(Near-field Acoustic Holography)算法的一种，类似声辐射模态（Acoustic Radiation Modes）法、快速傅里叶变换法（Fast Fourier Transform），边界元法（Boundary Element Method）等。可用于声场还原，噪声源识别定位，在振动与噪声控制、汽车、船舶、航空等领域等工程领域都有广泛应用等等。
基本理论代码部分
这部分主要用两种声源作为案例进行分析：单极子声源和简支板声源
ESM-主程序
tic
clear all
global k f
f_vec = 1000;%分析频率
% f_vec = 100:100:3000;
for nn = 1:length(f_vec)
f = f_vec(nn)% 激振频率
c = 343;
w = 2*pi*f;
k = w/c;
rho = 1.21;
F = 1;% 激振力
Lx = 0.5;
Ly = 0.5;
xy0 = [Lx/2,Ly/2];% 激振点位置
h = 0.001;%板厚度
z0 = 0;% 板子在z轴上的位置
PlateSize = [Lx,Ly];% 板子的尺寸
Ne = [121,121];% x轴、y轴方向上的板子划分成的单元格数量
%% 测量面
dxyh = 0.025;
dxyr = 0.025;
dxys = 0.025;
[XH,YH,ZH] = XYz_function(dxyh,0,Lx,0,Lx,0.1);
%XYz_function()是自定义的坐标位置网格函数，详见：函数一
[XR,YR,ZR] = XYz_function(dxyr,0,Lx,0,Lx,0.05);
[XS,YS,ZS] = XYz_function(dxys,0,Lx,0,Lx,-0.01);
Rs = [XS,YS,ZS];
Rr = [XR,YR,ZR];
Rh = [XH,YH,ZH];
%% 单极子(四个)
% rs11 = [-0.06 0 0];rs12 = [0.06 0 0];
% [Ph11,~] = ESM_point_source(Rh,rs11);
%ESM_point_source()单极子声源，详见：函数二
% [Ph12,~]= ESM_point_source(Rh,rs12);
% [Pth11,~]= ESM_point_source(Rr,rs11);
% [Pth12,~]= ESM_point_source(Rr,rs12);
% Prth = Pth11+Pth12;
% PH = Ph11+Ph12;
%% 简支板
[PH,~] = SimplySupportedPlate(f,F,xy0,PlateSize,h,Ne,Rh,z0);
%ESM_point_source()简支板声源，详见：函数三
[Prth,~] = SimplySupportedPlate(f,F,xy0,PlateSize,h,Ne,Rr,z0);
PH = awgn(PH,30,'measured');%添加高斯白噪声
%% 格林函数（传递函数）
Gpsh = CreatTM_p(Rh,Rs);
%CreatTM_p()格林函数/传递函数，详见：函数四
Gpsr = CreatTM_p(Rr,Rs);
[u,s,v] = csvd(Gpsh);
%csvd()奇异值分解（在这里处理对象为复数对，因此是csvd,而非svd）
lambda = gcv(u,s,PH);
%lambda = l_curve(u,s,PH);%L-曲线法
%gcv()广义交叉验证，与L-曲线法功能相同，仿真与实验各有所长，详见：函数五
Qreg = tikhonov(u,s,v,PH,lambda);
%tikhonov()基霍诺夫正则化，稳定求解过程，解决方程求逆时噪声放大问题，详见：函数六
% Qreg = Gpsh\PH;
Pr = Gpsr*Qreg;
%% 误差计算
Error_p(nn) = norm(abs(Pr)-abs(reshape(Prth,[],1)))...
/norm(abs(reshape(Prth,[],1)))*100
end
%% 误差绘图
figure(1)
plot(f_vec,Error_p,'ko-');
xlabel(' 频率(Hz)','Fontname','宋体','Fontsize',14);
ylabel('误差(%)','Fontname','宋体','Fontsize',14);
title('声压值误差曲线','Fontname','宋体','Fontsize',14);
box off;
%% 全息面（测量面）、重建面、重建面理论值声压值画图
figure(2)
subplot 131
contourf(abs(reshape(Prth,Lx/dxyh+1,Ly/dxyh+1)));
title('(a)理论值')
xlabel('x轴（测点数）');ylabel('y轴（测点数）');
hco = colorbar ;t = get(hco,'YTickLabel');t = strcat(t,'Pa');set(hco,'YTickLabel',t);
subplot 132
contourf(abs(reshape(PH,Lx/dxyr+1,Ly/dxyr+1)));
title('(b)测量值')
xlabel('x轴（测点数）');ylabel('y轴（测点数）');
hco = colorbar ;t = get(hco,'YTickLabel');t = strcat(t,'Pa');set(hco,'YTickLabel',t);
subplot 133
contourf(abs(reshape(Pr,Lx/dxyh+1,Ly/dxyh+1)));
title('(c)重建值')
xlabel('x轴（测点数）');ylabel('y轴（测点数）');
hco = colorbar ;t = get(hco,'YTickLabel');t = strcat(t,'Pa');set(hco,'YTickLabel',t);
函数一：绘制三维坐标系网格
说明：这里仅仅是笛卡尔坐标系下的三维坐标系，仅仅是研究的平面声全息，当然也可以把分析面设计成球体面、圆柱体面、锥形面等等，理论上，等效源法适用于任意形状面的声源所产生的声场。
function [X,Y,z]=XYz_function(dx,A,B,C,D,E)
x = A:dx:B;
y = C:dx:D;
[X,Y] = meshgrid(x,y);
X=reshape(X,[],1);
Y=reshape(Y,[],1);
z=reshape(E*ones(length(x),length(y)),[],1);
函数二：单极子声源声压值振速值
说明：在这部分，要说明单极子声源是理想的点声源，与脉动求声源的区别在于，单极子声源没有半径。
%%模拟单极子源所产生的声场（声压值和空气粒子振速值）
function [P,V] = ESM_point_source(Rm,rs)
%%基本参数
global f
c = 343;
rho = 1.21;
w = 2*pi*f;
k = w/c;
Q0 = 1*10^(-5);%仿真中单极子声源的强度系数
xs = rs(1);
ys = rs(2);
zs = rs(3);
Xh = Rm(:,1);
Yh = Rm(:,2);
Zh = Rm(:,3);
r1=sqrt((Xh-xs).^2+(Yh-ys).^2+(Zh-zs).^2);
P = (-1i*rho*w*Q0)./(4*pi*r1).*exp(1i*k*r1);%点声源声压公式频域
cosz = (Zh-zs)./r1;
V1 = (-1i*k*Q0)./(4*pi*r1);
V2 = 1-1./(1i*k*r1);
V3 = V1.*V2;
V4 = V3.*exp(1i*k*r1);
V = V4.*cosz; %点声源声压公式频域
% surf(abs(P));
函数三：简支板声源声压值、振速值
说明：简支板声源函数。包括声源的长宽厚尺寸、仿真种板子的划分单元数，板子密度（材料）、激励（激振器激励）位置、等都可以自己进行定义
function [p,v] = SimplySupportedPlate(f,F,xy0,PlateSize,h,Ne,Rm,z0)
%%%%% Simply Suported Aluminum Plate, point excited , origin is located at
%%%%% the left lower corner of plate.
%% excitation force
% f -- frequency of excitation force 激振频率
% F -- excitation force 激振力
% xy0 = （x,y,0) -- position of the excitation force 激振点位置
% h -- thickness of the plate 板厚度
% z0 -- 板子在z轴上的位置
% PlateSize = [size_x,size_y] -- 板子的尺寸
% Ne = [N_x,N_y] -- x轴、y轴方向上的板子划分成的单元格数量
% Rm = [XH;YH;ZH] -- 全息面位置向量
%% basic acoustic parameters
c = 343;
rho = 1.21;
w = 2*pi*f;
k = w/c;
%% basic parameters of the plate (Aluminum) 铝
% E = 7.1*10^10; %%%% Young's modulus 杨氏模量
% nju = 0.33; %%%% Possion's ratio 泊松比
% rhos = 2.7*10^3; %%%% density of the plate 板子密度
% D = E*h^3/(12*(1-nju^2)); %%% flexural rigidity of the plate 板的抗弯刚度
% alpha = (D/rhos/h)^(1/4);
%% basic parameters of the plate (Steel) 钢
E = 2.1*10^11; %%%% Young's modulus
nju = 0.28; %%%% Possion's ratio
rhos = 7.85*10^3; %%%% density of the plate
D = E*h^3/(12*(1-nju^2)); %%% flexural rigidity of the plate
alpha = (D/rhos/h)^(1/4);
%%
Lx = PlateSize(1);% x方向长度
Ly = PlateSize(2);% y方向长度
Nx = Ne(1);% x方向点数
Ny = Ne(2);% y方向点数
XH = Rm(:,1);% 测点x轴位置
YH = Rm(:,2);% 测点y轴位置
ZH = Rm(:,3);% 测点z轴位置
M = length(XH);
%% details of the measuring plane
if min(ZH-z0) > 0 %%% Rayleigh's Intergral
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
x = 0:dx:Lx;
y = 0:dy:Ly;
[X,Y] = meshgrid(x,y);
X = reshape(X,[],1);
Y = reshape(Y,[],1);
Wmn = zeros(length(x)*length(y),1);
for m = 1:20
for n = 1:20
PSImn0 = 2/sqrt(Lx*Ly)*sin(m*pi*xy0(1)/Lx).*sin(n*pi*xy0(2)/Ly);
PSImn = 2/sqrt(Lx*Ly)*sin(m*pi*X/Lx).*sin(n*pi*Y/Ly);
wmn = alpha^2*((m*pi/Lx)^2+(n*pi/Ly)^2);
Wmn = Wmn+(-F/(rhos*h))*PSImn0*PSImn/(w^2-wmn^2);
end
end
p = zeros(M,1);
v = zeros(M,1);
for ii = 1:M
r = sqrt((XH(ii)-X).^2+(YH(ii)-Y).^2+(ZH(ii)-z0)^2);
p(ii) = -w^2*rho/(2*pi)*sum(Wmn.*exp(1i*k*r)./r*dx*dy);
v(ii) = -1i*w/(2*pi)*sum(Wmn.*(1i*k*r-1)./r.^2.*exp(1i*k*r)*(ZH(ii)-z0)./r*dx*dy);
end
else %%% on the surface of the plate, only velocy is derived from displacement
Wmn = zeros(M,1);
for m = 1:20
for n = 1:20
PSImn0 = 2/sqrt(Lx*Ly)*sin(m*pi*xy0(1)/Lx).*sin(n*pi*xy0(2)/Ly);
PSImn = 2/sqrt(Lx*Ly)*sin(m*pi*XH/Lx).*sin(n*pi*YH/Ly);
wmn = alpha^2*((m*pi/Lx)^2+(n*pi/Ly)^2);
Wmn = Wmn+(-F/(rhos*h))*PSImn0*PSImn/(w^2-wmn^2);
end
end
v = -1i*w*Wmn;
p = NaN;
end
% p = reshape(p,length(xh),length(yh));
% v = reshape(v,length(xh),length(yh));
% pcolor(xh,yh,abs(p));;shading interp;
% figure(2); pcolor(xh,yh,abs(v));;shading interp;
函数四：传递函数
说明：这部分是针对传递函数这里也称格林函数，进行代码实现部分，没啥好说的，应该不难理解
function Gp = CreatTM_p(Rh,Rs)
% Rh为全息面位置向量，共三列，每行 【x,y,z】
global k
Xh = Rh(:,1);
Yh = Rh(:,2);
Zh = Rh(:,3);
Xs = Rs(:,1);
Ys = Rs(:,2);
Zs = Rs(:,3);
[XS,XH] = meshgrid(Xs,Xh);
[YS,YH] = meshgrid(Ys,Yh);
[ZS,ZH] = meshgrid(Zs,Zh);
R_hs = sqrt((XH-XS).^2+(YH-YS).^2+(ZH-ZS).^2);
Gp = exp(1i*k*R_hs)./(4*pi*R_hs);
函数五：正则化参数选取
说明：GCV广义交叉验证法，功能：在对方程进行正则化求解过程中，选取正则化参数的一种方法，与L-曲线法类似，但是GCV方法有时候更适用于实验室实验的计算，而L-曲线法在数值仿真中的计算效果更好。关于为什么要使用正则化这个问题，简单地说是为了抑制噪声，详细的请自行参阅相关文献，资源很多，不再累述。
另外：这段函数还有另外一个小函数模块，会在这段代码后面补上。
function [reg_min,G,reg_param] = gcv(U,s,b,method)
%GCV Plot the GCV function and find its minimum.
%
% [reg_min,G,reg_param] = gcv(U,s,b,method)
% [reg_min,G,reg_param] = gcv(U,sm,b,method) , sm = [sigma,mu]
%
% Plots the GCV-function
% || A*x - b ||^2
% G = -------------------
% (trace(I - A*A_I)^2
% as a function of the regularization parameter reg_param.
% Here, A_I is a matrix which produces the regularized solution.
%
% The following methods are allowed:
% method = 'Tikh' : Tikhonov regularization (solid line )
% method = 'tsvd' : truncated SVD or GSVD (o markers )
% method = 'dsvd' : damped SVD or GSVD (dotted line)
% If method is not specified, 'Tikh' is default.
%
% If any output arguments are specified, then the minimum of G is
% identified and the corresponding reg. parameter reg_min is returned.
% Per Christian Hansen, IMM, Dec. 16, 2003.
% Reference: G. Wahba, "Spline Models for Observational Data",
% SIAM, 1990.
% Set defaults.
if (nargin==3), method='Tikh'; end % Default method.
npoints = 200; % Number of points on the curve.
smin_ratio = 16*eps; % Smallest regularization parameter.
% Initialization.
[m,n] = size(U); [p,ps] = size(s);
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==2)
s = s(p:-1:1,1)./s(p:-1:1,2); beta = beta(p:-1:1);
end
if (nargout > 0), find_min = 1; else find_min = 0; end
if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))
% Vector of regularization parameters.
reg_param = zeros(npoints,1); G = reg_param; s2 = s.^2;
reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
% Intrinsic residual.
delta0 = 0;
if (m > n & beta2 > 0), delta0 = beta2; end
% Vector of GCV-function values.
for i=1:npoints
G(i) = gcvfun(reg_param(i),s2,beta(1:p),delta0,m-n);
end
% Plot GCV function.
loglog(reg_param,G,'-'), xlabel('\lambda'), ylabel('G(\lambda)')
title('GCV function')
% Find minimum, if requested.
if (find_min)
[minG,minGi] = min(G); % Initial guess.
reg_min = fminbnd('gcvfun',...
reg_param(min(minGi+1,npoints)),reg_param(max(minGi-1,1)),...
optimset('Display','off'),s2,beta(1:p),delta0,m-n); % Minimizer.
minG = gcvfun(reg_min,s2,beta(1:p),delta0,m-n); % Minimum of GCV function.
ax = axis;
HoldState = ishold; hold on;
loglog(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
title(['GCV function, minimum at \lambda = ',num2str(reg_min)])
axis(ax)
if (~HoldState), hold off; end
end
elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))
% Vector of GCV-function values.
rho2(p-1) = abs(beta(p))^2;
if (m > n & beta2 > 0), rho2(p-1) = rho2(p-1) + beta2; end
for k=p-2:-1:1, rho2(k) = rho2(k+1) + abs(beta(k+1))^2; end
for k=1:p-1
G(k) = rho2(k)/(m - k + (n - p))^2;
end
reg_param = [1:p-1]';
% Plot GCV function.
semilogy(reg_param,G,'o'), xlabel('k'), ylabel('G(k)')
title('GCV function')
% Find minimum, if requested.
if (find_min)
[minG,reg_min] = min(G);
ax = axis;
HoldState = ishold; hold on;
semilogy(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
title(['GCV function, minimum at k = ',num2str(reg_min)])
axis(ax);
if (~HoldState), hold off; end
end
elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))
% Vector of regularization parameters.
reg_param = zeros(npoints,1); G = reg_param;
reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
% Intrinsic residual.
delta0 = 0;
if (m > n & beta2 > 0), delta0 = beta2; end
% Vector of GCV-function values.
for i=1:npoints
G(i) = gcvfun(reg_param(i),s,beta(1:p),delta0,m-n,1);
end
% Plot GCV function.
loglog(reg_param,G,':'), xlabel('\lambda'), ylabel('G(\lambda)')
title('GCV function')
% Find minimum, if requested.
if (find_min)
[minG,minGi] = min(G); % Initial guess.
reg_min = fminbnd('gcvfun',...
reg_param(min(minGi+1,npoints)),reg_param(max(minGi-1,1)),...
optimset('Display','off'),s,beta(1:p),delta0,m-n,1); % Minimizer.
minG = gcvfun(reg_min,s,beta(1:p),delta0,m-n,1); % Minimum of GCV function.
ax = axis;
HoldState = ishold; hold on;
loglog(reg_min,minG,'*r',[reg_min,reg_min],[minG/1000,minG],':r')
title(['GCV function, minimum at \lambda = ',num2str(reg_min)])
axis(ax)
if (~HoldState), hold off; end
end
elseif (strncmp(method,'mtsv',4) | strncmp(method,'ttls',4))
error('The MTSVD and TTLS methods are not supported')
else, error('Illegal method'), end
gcvfun
说明：函数五：gcv的内置函数
function G = gcvfun(lambda,s2,beta,delta0,mn,dsvd)
% Auxiliary routine for gcv. PCH, IMM, 12/29/97.
if (nargin==5)
f = (lambda^2)./(s2 + lambda^2);
else
f = lambda./(s2 + lambda);
end
G = (norm(f.*beta)^2 + delta0)/(mn + sum(f))^2;
L-曲线法：选取正则化参数
说明：在正则化参数选取过程中的另外一总算法，这些，我只是拿过来用，里面的很多东西我也不懂，其程序以及内置函数如下：
function [reg_corner,rho,eta,reg_param] = l_curve(U,sm,b,method,L,V)
%L_CURVE Plot the L-curve and find its "corner".
%
% [reg_corner,rho,eta,reg_param] =
% l_curve(U,s,b,method)
% l_curve(U,sm,b,method) , sm = [sigma,mu]
% l_curve(U,s,b,method,L,V)
%
% Plots the L-shaped curve of eta, the solution norm || x || or
% semi-norm || L x ||, as a function of rho, the residual norm
% || A x - b ||, for the following methods:
% method = 'Tikh' : Tikhonov regularization (solid line )
% method = 'tsvd' : truncated SVD or GSVD (o markers )
% method = 'dsvd' : damped SVD or GSVD (dotted line)
% method = 'mtsvd' : modified TSVD (x markers )
% The corresponding reg. parameters are returned in reg_param. If no
% method is specified then 'Tikh' is default. For other methods use plot_lc.
%
% Note that 'Tikh', 'tsvd' and 'dsvd' require either U and s (standard-
% form regularization) or U and sm (general-form regularization), while
% 'mtvsd' requires U and s as well as L and V.
%
% If any output arguments are specified, then the corner of the L-curve
% is identified and the corresponding reg. parameter reg_corner is
% returned. Use routine l_corner if an upper bound on eta is required.
%
% If the Spline Toolbox is not available and reg_corner is requested,
% then the routine returns reg_corner = NaN for 'tsvd' and 'mtsvd'.
% Reference: P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
% the regularization of discrete ill-posed problems", SIAM J. Sci.
% Comput. 14 (1993), pp. 1487-1503.
% Per Christian Hansen, IMM, Sep. 13, 2001.
% Set defaults.
if (nargin==3), method='Tikh'; end % Tikhonov reg. is default.
npoints = 200; % Number of points on the L-curve for Tikh and dsvd.
smin_ratio = 16*eps; % Smallest regularization parameter.
% Initialization.
[m,n] = size(U); [p,ps] = size(sm);
if (nargout > 0), locate = 1; else locate = 0; end
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==1)
s = sm; beta = beta(1:p);
else
s = sm(p:-1:1,1)./sm(p:-1:1,2); beta = beta(p:-1:1);
end
xi = beta(1:p)./s;
if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))
eta = zeros(npoints,1); rho = eta; reg_param = eta; s2 = s.^2;
reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
for i=1:npoints
f = s2./(s2 + reg_param(i)^2);
eta(i) = norm(f.*xi);
rho(i) = norm((1-f).*beta(1:p));
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
marker = '-'; pos = .8; txt = 'Tikh.';
elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))
eta = zeros(p,1); rho = eta;
eta(1) = abs(xi(1))^2;
for k=2:p, eta(k) = eta(k-1) + abs(xi(k))^2; end
eta = sqrt(eta);
if (m > n)
if (beta2 > 0), rho(p) = beta2; else rho(p) = eps^2; end
else
rho(p) = eps^2;
end
for k=p-1:-1:1, rho(k) = rho(k+1) + abs(beta(k+1))^2; end
rho = sqrt(rho);
reg_param = [1:p]'; marker = 'o'; pos = .75;
if (ps==1)
U = U(:,1:p); txt = 'TSVD';
else
U = U(:,1:p); txt = 'TGSVD';
end
elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))
eta = zeros(npoints,1); rho = eta; reg_param = eta;
reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
for i=1:npoints
f = s./(s + reg_param(i));
eta(i) = norm(f.*xi);
rho(i) = norm((1-f).*beta(1:p));
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
marker = ':'; pos = .85;
if (ps==1), txt = 'DSVD'; else txt = 'DGSVD'; end
elseif (strncmp(method,'mtsv',4))
if (nargin~=6)
error('The matrices L and V must also be specified')
end
[p,n] = size(L); rho = zeros(p,1); eta = rho;
[Q,R] = qr(L*V(:,n:-1:n-p),0);
for i=1:p
k = n-p+i;
Lxk = L*V(:,1:k)*xi(1:k);
zk = R(1:n-k,1:n-k)\(Q(:,1:n-k)'*Lxk); zk = zk(n-k:-1:1);
eta(i) = norm(Q(:,n-k+1:p)'*Lxk);
if (i < p)
rho(i) = norm(beta(k+1:n) + s(k+1:n).*zk);
else
rho(i) = eps;
end
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
reg_param = [n-p+1:n]'; txt = 'MTSVD';
U = U(:,reg_param); sm = sm(reg_param);
marker = 'x'; pos = .7; ps = 2; % General form regularization.
else, error('Illegal method'), end
% Locate the "corner" of the L-curve, if required. If the Spline
% Toolbox is not available, return NaN for reg_corner.
if (locate)
SkipCorner = ( (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4) | ...
strncmp(method,'mtsv',4)) & exist('splines')~=7 );
if (SkipCorner)
reg_corner = NaN;
else
[reg_corner,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,sm,b,method);
end
end
% Make plot.
plot_lc(rho,eta,marker,ps,reg_param);
if (locate & ~SkipCorner)
ax = axis;
HoldState = ishold; hold on;
loglog([min(rho)/100,rho_c],[eta_c,eta_c],':r',...
[rho_c,rho_c],[min(eta)/100,eta_c],':r')
title(['L-curve, ',txt,' corner at ',num2str(reg_corner)]);
axis(ax)
if (~HoldState), hold off; end
end
l_corner
function [reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M)
%L_CORNER Locate the "corner" of the L-curve.
%
% [reg_c,rho_c,eta_c] =
% l_corner(rho,eta,reg_param)
% l_corner(rho,eta,reg_param,U,s,b,method,M)
% l_corner(rho,eta,reg_param,U,sm,b,method,M) , sm = [sigma,mu]
%
% Locates the "corner" of the L-curve in log-log scale.
%
% It is assumed that corresponding values of || A x - b ||, || L x ||,
% and the regularization parameter are stored in the arrays rho, eta,
% and reg_param, respectively (such as the output from routine l_curve).
%
% If nargin = 3, then no particular method is assumed, and if
% nargin = 2 then it is issumed that reg_param = 1:length(rho).
%
% If nargin >= 6, then the following methods are allowed:
% method = 'Tikh' : Tikhonov regularization
% method = 'tsvd' : truncated SVD or GSVD
% method = 'dsvd' : damped SVD or GSVD
% method = 'mtsvd' : modified TSVD,
% and if no method is specified, 'Tikh' is default. If the Spline Toolbox
% is not available, then only 'Tikh' and 'dsvd' can be used.
%
% An eighth argument M specifies an upper bound for eta, below which
% the corner should be found.
% The following functions from the Spline Toolbox are needed if
% method differs from 'Tikh' or 'dsvd':
% fnder, ppbrk, ppmak, ppual, sp2pp, sorted, spbrk, spmak, sprpp.
% Per Christian Hansen, IMM, Dec. 12, 2002.
% Set default regularization method.
if (nargin <= 3)
method = 'none';
if (nargin==2), reg_param = [1:length(rho)]'; end
else
if (nargin==6), method = 'Tikh'; end
end
% Set threshold for skipping very small singular values in the
% analysis of a discrete L-curve.
s_thr = eps; % Neglect singular values less than s_thr.
% Set default parameters for treatment of discrete L-curve.
deg = 2; % Degree of local smooting polynomial.
q = 2; % Half-width of local smoothing interval. orig = 2
order = 4; % Order of fitting 2-D spline curve.
% Initialization.
if (length(rho) < order)
error('Too few data points for L-curve analysis')
end
if (nargin > 3)
[p,ps] = size(s); [m,n] = size(U);
beta = U'*b;
if (m>n), b0 = b - U*beta; end
if (ps==2)
s = s(p:-1:1,1)./s(p:-1:1,2);
U = U(:,p:-1:1);
beta = beta(p:-1:1);
end
xi = beta./s;
end
% Restrict the analysis of the L-curve according to M (if specified).
if (nargin==8)
index = find(eta < M);
rho = rho(index); eta = eta(index); reg_param = reg_param(index);
end
if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))
% The L-curve is differentiable; computation of curvature in
% log-log scale is easy.
% Compute g = - curvature of L-curve.
g = lcfun(reg_param,s,beta,xi);
% Locate the corner. If the curvature is negative everywhere,
% then define the leftmost point of the L-curve as the corner.
[gmin,gi] = min(g);
reg_c = fminbnd('lcfun',...
reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
optimset('Display','off'),s,beta,xi); % Minimizer.
kappa_max = - lcfun(reg_c,s,beta,xi); % Maximum curvature.
if (kappa_max < 0)
lr = length(rho);
reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
else
f = (s.^2)./(s.^2 + reg_c^2);
eta_c = norm(f.*xi);
rho_c = norm((1-f).*beta);
if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
end
elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4) | ...
strncmp(method,'mtsv',4) | strncmp(method,'none',4))
% The L-curve is discrete and may include unwanted fine-grained
% corners. Use local smoothing, followed by fitting a 2-D spline
% curve to the smoothed discrete L-curve.
% Check if the Spline Toolbox exists, otherwise return.
if (exist('splines')~=7)
error('The Spline Toolbox in not available so l_corner cannot be used')
end
% For TSVD, TGSVD, and MTSVD, restrict the analysis of the L-curve
% according to s_thr.
if (nargin > 3)
if (nargin==8) % In case the bound M is in action.
s = s(index,:);
beta = beta(index);
xi = xi(index);
end
index = find(s > s_thr);
rho = rho(index); eta = eta(index); reg_param = reg_param(index);
s = s(index); beta = beta(index); xi = xi(index);
end
% Convert to logarithms.
lr = length(rho);
lrho = log(rho); leta = log(eta); slrho = lrho; sleta = leta;
% For all interior points k = q+1:length(rho)-q-1 on the discrete
% L-curve, perform local smoothing with a polynomial of degree deg
% to the points k-q:k+q.
v = [-q:q]'; A = zeros(2*q+1,deg+1); A(:,1) = ones(length(v),1);
for j = 2:deg+1, A(:,j) = A(:,j-1).*v; end
for k = q+1:lr-q-1
cr = A\lrho(k+v); slrho(k) = cr(1);
ce = A\leta(k+v); sleta(k) = ce(1);
end
% Fit a 2-D spline curve to the smoothed discrete L-curve.
sp = spmak([1:lr+order],[slrho';sleta']);
pp = ppbrk(sp2pp(sp),[4,lr+1]);
% Extract abscissa and ordinate splines and differentiate them.
% Compute as many function values as default in spleval.
P = spleval(pp); dpp = fnder(pp);
D = spleval(dpp); ddpp = fnder(pp,2);
DD = spleval(ddpp);
ppx = P(1,:); ppy = P(2,:);
dppx = D(1,:); dppy = D(2,:);
ddppx = DD(1,:); ddppy = DD(2,:);
% Compute the corner of the discretized .spline curve via max. curvature.
% No need to refine this corner, since the final regularization
% parameter is discrete anyway.
% Define curvature = 0 where both dppx and dppy are zero.
k1 = dppx.*ddppy - ddppx.*dppy;
k2 = (dppx.^2 + dppy.^2).^(1.5);
I_nz = find(k2 ~= 0);
kappa = zeros(1,length(dppx));
kappa(I_nz) = -k1(I_nz)./k2(I_nz);
[kmax,ikmax] = max(kappa);
x_corner = ppx(ikmax); y_corner = ppy(ikmax);
% Locate the point on the discrete L-curve which is closest to the
% corner of the spline curve. Prefer a point below and to the
% left of the corner. If the curvature is negative everywhere,
% then define the leftmost point of the L-curve as the corner.
if (kmax < 0)
reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
else
index = find(lrho < x_corner & leta < y_corner);
if (length(index) > 0)
[dummy,rpi] = min((lrho(index)-x_corner).^2 + (leta(index)-y_corner).^2);
rpi = index(rpi);
else
[dummy,rpi] = min((lrho-x_corner).^2 + (leta-y_corner).^2);
end
reg_c = reg_param(rpi); rho_c = rho(rpi); eta_c = eta(rpi);
end
elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))
% The L-curve is differentiable; computation of curvature in
% log-log scale is easy.
% Compute g = - curvature of L-curve.
g = lcfun(reg_param,s,beta,xi,1);
% Locate the corner. If the curvature is negative everywhere,
% then define the leftmost point of the L-curve as the corner.
[gmin,gi] = min(g);
reg_c = fminbnd('lcfun',...
reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
optimset('Display','off'),s,beta,xi,1); % Minimizer.
kappa_max = - lcfun(reg_c,s,beta,xi,1); % Maximum curvature.
if (kappa_max < 0)
lr = length(rho);
reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
else
f = s./(s + reg_c);
eta_c = norm(f.*xi);
rho_c = norm((1-f).*beta);
if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
end
else, error('Illegal method'), end
lsolve
function x = lsolve(L,y,W,T)
%LSOLVE Utility routine for "preconditioned" iterative methods.
% x = lsolve(L,y,W,T)
% Computes the vector
% x = L_p*y
% where L_p is the A-weighted generalized inverse of L.
% Here, L is a p-by-n matrix, W holds a basis for the null space of L,
% and T is a utility matrix which should be computed by routine pinit.
% Alternatively, L is square, and W and T are not needed.
% Notice that x and y may be matrices, in which case
% x(:,i) = L_p*y(:,i) .
% Reference: P. C. Hansen, "Rank-Deficient and Discrete Ill-Posed Problems.
% Numerical Aspects of Linear Inversion", SIAM, Philadelphia, 1997.
% Per Christian Hansen, IMM, 07/29/97.
% Initialization.
[p,n] = size(L); nu = n-p; ly = size(y,2);
% Special treatment of square L.
if (nu==0), x = L\y; return; end
% The general case.
x = L(:,1:p)\y;
x = [x;zeros(nu,ly)] - W*(T(:,1:p)*x);
函数六：正则化稳定求解过曾
说明：基霍诺夫正则化
function [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)
%TIKHONOV Tikhonov regularization.
%
% [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)
% [x_lambda,rho,eta] = tikhonov(U,sm,X,b,lambda,x_0) , sm = [sigma,mu]
%
% Computes the Tikhonov regularized solution x_lambda. If the SVD
% is used, i.e. if U, s, and V are specified, then standard-form
% regularization is applied:
% min { || A x - b ||^2 + lambda^2 || x - x_0 ||^2 } .
% If, on the other hand, the GSVD is used, i.e. if U, sm, and X are
% specified, then general-form regularization is applied:
% min { || A x - b ||^2 + lambda^2 || L (x - x_0) ||^2 } .
%
% If x_0 is not specified, then x_0 = 0 is used.
%
% Note that x_0 cannot be used if A is underdetermined and L ~= I.
%
% If lambda is a vector, then x_lambda is a matrix such that
% x_lambda = [ x_lambda(1), x_lambda(2), ... ] .
%
% The solution norm(standard-form case) or seminorm (general-form
% case) and the residual norm are returned in eta and rho.
% Per Christian Hansen, IMM, April 14, 2003.
% Reference: A. N. Tikhonov & V. Y. Arsenin, "Solutions of
% Ill-Posed Problems", Wiley, 1977.
% Initialization.
if (min(lambda)<0)
error('Illegal regularization parameter lambda')
end
m = size(U,1);
n = size(V,1);
[p,ps] = size(s);
beta = U(:,1:p)'*b;
zeta = s(:,1).*beta;
ll = length(lambda); x_lambda = zeros(n,ll);
rho = zeros(ll,1); eta = zeros(ll,1);
% Treat each lambda separately.
if (ps==1)
% The standard-form case.
if (nargin==6), omega = V'*x_0; end
for i=1:ll
if (nargin==5)
x_lambda(:,i) = V(:,1:p)*(zeta./(s.^2 + lambda(i)^2));
rho(i) = lambda(i)^2*norm(beta./(s.^2 + lambda(i)^2));
else
x_lambda(:,i) = V(:,1:p)*...
((zeta + lambda(i)^2*omega)./(s.^2 + lambda(i)^2));
rho(i) = lambda(i)^2*norm((beta - s.*omega)./(s.^2 + lambda(i)^2));
end
eta(i) = norm(x_lambda(:,i));
end
if (nargout > 1 & size(U,1) > p)
rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);
end
elseif (m>=n)
% The overdetermined or square general-form case.
gamma2 = (s(:,1)./s(:,2)).^2;
if (nargin==6), omega = V\x_0; omega = omega(1:p); end
if (p==n)
x0 = zeros(n,1);
else
x0 = V(:,p+1:n)*U(:,p+1:n)'*b;
end
for i=1:ll
if (nargin==5)
xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);
x_lambda(:,i) = V(:,1:p)*xi + x0;
rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));
else
xi = (zeta + lambda(i)^2*(s(:,2).^2).*omega)./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);
x_lambda(:,i) = V(:,1:p)*xi + x0;
rho(i) = lambda(i)^2*norm((beta - s(:,1).*omega)./(gamma2 + lambda(i)^2));
end
eta(i) = norm(s(:,2).*xi);
end
if (nargout > 1 & size(U,1) > p)
rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);
end
else
% The underdetermined general-form case.
gamma2 = (s(:,1)./s(:,2)).^2;
if (nargin==6), error('x_0 not allowed'), end
if (p==m)
x0 = zeros(n,1);
else
x0 = V(:,p+1:m)*U(:,p+1:m)'*b;
end
for i=1:ll
xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);
x_lambda(:,i) = V(:,1:p)*xi + x0;
rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));
eta(i) = norm(s(:,2).*xi);
end
end
声名
1.在学习的过程中，往往都是通过查阅论文进行研究，而我发现在这个领域的研究成果一般都是闭元的，很少有具体的实验案例，毕竟这个学科领域正处在理论研究阶段，大多数的研究案例都是通过数学公式将理论进行阐释，对于很多初学者来讲，往往最容易带来困惑的，往往都是一些细节。而在互联网上的一些贴吧、论坛、网站进行学术交流时，效率实在是慢，各大门户的网站各具特色，同时也带来信息共享渠道不兼容等问题，总之吧，我写这篇博客纯属学术交流，为初学者提供一点帮助。
2.编程仿真都是通过Matlab环境进行编写，实验验证也是通过Matlab进行的，在这里分享的程序也都是基于Matlab环境下的代码。
3.请转载时注明来源。