%本程序可以實現用svd分解法和阻尼最小二乘法來進行反演，如果choose=0，實現svd分解法，如果choose=1，實現阻尼最小二乘法
function [G,D,DD,M,MM]=fanyan
clear;
x=0:0.02:1;
y=0:1:10;
n=length(x);
m=length(y);
G=zeros(m,n);
choose=input('if you want to use SVD method,choose=0,then damped least-square method,choose=1,so choose=?')
zuni=0.0001;
pi=3.14;
%計算核函數
for i=1:m
    for j=1:n
        G(i,j)=exp(-i*x(j));
    end
end
%給定理論模型
M=zeros(n,1);
for i=1:n
    M(i,1)=1-0.5*cos(2*pi*x(i));
end
D=G*M;
%在理論計算的數據上加噪聲得到觀測值
N=0.001*randn(m,1);%因為噪聲具有隨機性，所以得出的結果也可能不一樣
DD=D+N; 
%兩種反演方法的選擇
[U,S,V]=svd(G);
 if   choose==0
     S(11,11)=0; 
     S(10,10)=0;
     S(9,9)=0;
     S(8,8)=0;
     S(7,7)=0;
     SS=pinv(S);
 else
     [U1,S1,V1]=svd(G');
    SS=pinv(S1*S+zuni*eye(n,n))*S1;
 end   
MM=V*SS*U'*DD;%計算反演模型
figure;
plot(x,M,x,MM,'r-*');
legend('理論模型曲線','反演模型曲線');