此函数用FISTA算法解决压缩感知

此函数用FISTA算法解决压缩感知，注释很详细，里面包含测试图片

%  此函数用FISTA算法解决压缩感知（傅里叶变换做测量，磁共振应用）
%  参考文献：Amir Beck and Marc Teboulle “A Fast Iterative Shrinkage-Thresholding Algorithm for 
%  Linear Inverse Problems“ SIAM J. IMAGING SCIENCES Vol. 2 No. 1 pp.183-202 2009.
%  编程人：沙威 香港大学
%  Email: wsha@eee.hku.hk; dr.weisha@gmail.com
%  转载时请保留上面注解

function FISTA_FT;

clc;clear

%  图像库
A={‘phantom256.bmp‘ ‘fruits256.bmp‘ ‘cameraman256.bmp‘ ‘lena256.bmp‘...
    ‘peppers256.bmp‘ ‘boat256.bmp‘ ‘baboon256.bmp‘};

I=imread（A{3}）;           %  1表示最经典的phantom图像，可以尝试其它的图像
X=double（I）/255;          %  原始矩阵归一化
load ww                   %  小波变换矩阵
W_i=WL_for（wwX）;         %  小波变换
max_e=max（max（full（abs（W_i））））; %  最大系数
threshold=0.005;          %  稀疏化阈值（小波变换）
th_v=（max_e*threshold）;   %  阈值

N=length（X）;              %  原始图像大小

% 空域随机矩阵（Y=A*X测量，A病态，需要用共轭梯度技术，需要修正）
% A=rand（N/4N）;          %  随机测量矩阵（25%测量）
% A=A/norm（A）;            %  测量矩阵归一化

load mask_radial;         %  mask矩阵（30%测量）
Y=FT_for（mask_matrixX）;  %  傅里叶变换系数作为测量

X_r=zeros（NN）;  %  恢复矩阵
N_iter=100;      %  迭代次数
c=1;             %  代理函数参数
gamma=1/c;       %  迭代系数（可以优化）

t=1;             %  自适应迭代系数

for m=1:N_iter
    N_iter-m
    
    t_p=t;                            %  记录之前迭代系数
    X_p=X_r;                          %  记录之前的小波域系数

    X_t=WL_back（wwX_r）;              %  小波反变换，得到时域图像
    X_m=FT_for（mask_matrixX_t）;      %  傅里叶变换域测量
    R_m=FT_back（mask_matrix（Y-X_m））; %  测量残差
    X_s=X_r+gamma*WL_for（wwR_m）;     %  小波正变换，更新小波域系数
    X_r=Threshold_F（X_sth_v/c）;      %  小波域系数软阈值

    t=（1+sqrt（1+4*t^2））/2;            %  新的迭代系数
    X_r=X_r+（t_p-1）/t*（X_r-X_p）;      %  更新小波域系数方法1
    %X_r=X_r+（t-1）/（t+1）*（X_r-X_p）;   %  更新小波域系数方法2
    
    th_v=th_v*0.9;                    %  软阈值松弛
    
    if （th_v<1e-3）                    %  阈值很小时候截断迭代
        break;
    end
    
    %  和原始图像相比的峰值信噪比PSNR
    errorx=sum（sum（abs（WL_back（wwX_r）-X）.^2））;   %  MSE误差
    psnr=10*log10（1*1/（errorx/N/N））;  %  PSNR
    disp（‘峰值信噪比：‘）
    disp（psnr）
 
end

figure（1）
subplot（221）
imshow（uint8（X*255））
title（‘原始图像‘）
subplot（222）
imshow（uint8（real（Y）*255））
title（‘傅里叶域测量‘）
subplot（223）
imshow（uint8（full（WL_back（wwX_r））*255））
ss=sprintf（‘恢复图像（PSNR: %0.2f）‘psnr）;
title（ss）
subplot（224）
imshow（uint8（full（WL_back（wwX_r）-X）*255））
title（‘误差图像‘）

%  小波变换和反变换
function MM=WL_for（wwM）
MM=（ww*（M）*ww‘）;     %  2d image
%MM=（ww*sparse（M））;  %  1d signal

function MM=WL_back（wwM）
MM=（ww‘*（M）*ww）;      %  2d image
%MM=（ww‘*sparse（M））;  %  1d signal

%  阈值函数
function X_th=Threshold_F（Xthreshold）

%  软阈值
flag1=X>（threshold）;
flag2=X<-（threshold）;
flag3=flag1 | flag2;
X_th=X-double（flag1）.*threshold;
X_th=X_th+double（flag2）.*threshold;
X_th=X_th.*flag3;

% 傅里叶变换和反变换
function MM=FT_for（maskM）
M=real（M）;
MM=mask.*（fftshift（fft2（M）））;

function MM=FT_back（maskM）
MM=real（ifft2（ifftshift（mask.*M）））;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      66614  2015-06-17 20:13  FISTA\baboon256.bmp

     文件      66614  2015-06-17 20:13  FISTA\boat256.bmp

     文件      66358  2015-06-17 20:13  FISTA\brain256.bmp

     文件      66614  2015-06-17 20:13  FISTA\cameraman256.bmp

     文件       3298  2015-06-17 20:13  FISTA\FISTA_FT.m

     文件       2907  2018-07-07 16:37  FISTA\FISTA_Radon.m

     文件      66614  2015-06-17 20:13  FISTA\fruits256.bmp

     文件      66614  2015-06-17 20:13  FISTA\lena256.bmp

     文件       6031  2015-06-17 20:13  FISTA\mask_radial.mat

     文件      66614  2015-06-17 20:13  FISTA\peppers256.bmp

     文件       3824  2015-06-17 20:13  FISTA\phantom256.bmp

     文件      80168  2015-06-17 20:13  FISTA\ww.mat

     目录          0  2018-07-06 15:13  FISTA

----------- ---------  ---------- -----  ----

               562270                    13