clear
clc
flag=2;  % flag=1, 向量; flag=2, 单变量(一个数值); flag=3, MATLAB自带的软阈值函数;  flag=4, 自定义软阈值函数
%% 向量的情况 2D shrinkage formula
% 最小化问题：z(x)=alpha*||x||+(beta/2)*||x-y||^2
% alpha=0.5;
% beta=1;
% 用的2范数
if flag==1
    % 法1：解析解
    % 当||y||>alpha，x=(||y||-alpha/beta).*(y/||y||)
    y=0.1:0.1:1;
    num=length(y);
    yy=norm(y, 2);
    x_1=((yy-0.5)/yy).*y;
    z_1=0.5*norm(x_1, 2)+0.5.*(x_1-y)*(x_1-y)';
    % 法2：数值解--只要重复次数够大，能逼近解析解的结果
    % x=(||y||/((alpha/beta)+||y||)).*y
    N=50; %重复迭代N次
    x_2=rand(1, num);
    for i=1:N
        xx=norm(x_2, 2);
        x_2=(xx/(0.5+xx)).*y;
    end
    z_2=0.5*norm(x_2, 2)+0.5.*(x_2-y)*(x_2-y)';
    % 结果
    x_result=[x_1; x_2] % x的结果, 第一行是方法1的结果, 第二行是方法2的结果
    z=[z_1; z_2] % 最小化问题z的结果, 第一行是方法1的结果, 第二行是方法2的结果
    
%% 单变量的情况 1D shrinkage formula
elseif flag==2
    % z(x)=(beta/2)*(x-y)^2+alpha*|x|
    % alpha=0.5, beta=1
    % 当|y|>alpha, x=(|y|-alpha)*(y/|y|)
    x=0.1:0.1:1;
    y=1;
    z=0.5.*(x-y).^2+0.5.*abs(x);
    [z_min, index]=min(z);
    plot(x, z);
    hold on
    plot(x(index), z_min, 'ro');
    xlabel('x');
    ylabel('z');
    x_result=x(index)
    z=z_min
%% MATLAB自带的软阈值函数  1D shrinkage formula
elseif flag==3
    x=-5:0.01:5;
    thr=2;
    %  z(x)=(beta/2)*(x-y)^2+alpha*|x|
    % thr=alpha/beta
    ysoft=wthresh(x,'s',thr); % 软阈值
    ythard=wthresh(x,'h',thr); % 硬阈值
    plot(x, x, 'b--', x, ythard, 'g-', x, ysoft, 'r-', 'LineWidth',1.3);
    legend('Original Signal','Hard Threshold','Soft Threshold', 'Location', 'northwest'); %图例的设置
    % 设置网格线
    grid on;
    set(gca, 'FontName','Times New Roman');
    set(gca, 'GridLineStyle', ':');
    set(gca, 'GridAlpha', 1);
    saveas(gcf,sprintf('matlab_shrinkage.jpg'),'bmp'); %保存图片
%% 自定义软阈值函数  1D shrinkage formula
elseif flag==4
    x=-5:0.01:5;
    thr=2;
    y=Soft_Threshold(x,thr);
    plot(x, y, 'r-', 'LineWidth',1.3)
    title('Soft Threshold');  
    % 设置网格线
    grid on;
    set(gca, 'FontName','Times New Roman');
    set(gca, 'GridLineStyle', ':');
    set(gca, 'GridAlpha', 1);
    saveas(gcf,sprintf('soft_shrinkage.jpg'),'bmp'); %保存图片
end
% 自定义软阈值函数
%  z(x)=(beta/2)*(x-y)^2+alpha*|x|
% thr=alpha/beta
function y=Soft_Threshold(x,thr)
    y=sign(x).*max(abs(x) - thr,0);
end