求用MP算法分解重建图像的matlab代码

这里是一维信号的重建：

%基于MP算法

clcclear

%观测向手野量y的长度M=80，即采样率M/N=0.3

N=256

K=15%信号稀疏度为15

M=80%

x = zeros(N,1)

q = randperm(N)

x(q(1:K)) =randn(K,1) %原始信号

%构造高斯测量矩阵，用以随机采样

Phi = randn(M,N)*sqrt(1/M)

for i = 1:N

Phi(:,i) = Phi(:,i)/norm(Phi(:,i))

end

y=Phi*x %获得线性测量

%用MP算法开始迭代重构

m=2*K %总的迭代次数

r_n=y % 残差值初始值

x_find=zeros(N,1) %x_find为MP算法恢复世罩的信号

for times=1:m

for col=1:N

neiji(col)=Phi(:,col)'*r_n %计算当前残差和感知矩阵每一列的内积

end

[val,pos]=max(abs(neiji)) %找出内积中绝对值最大的元素和它的对应的感知矩毕返喊阵的列pos

x_find(pos)=x_find(pos)+neiji(pos) %计算新的近似x_find

r_n=r_n-neiji(pos)*Phi(:,pos) %更新残差

end

subplot(3,1,1)plot(x)title('target')

subplot(3,1,2)plot(x_find)title('reconstruct')

subplot(3,1,3)plot(r_n)title('残差')

%用phantom函数可以获得仿体图像；

%用randon可获得不同角度的一维投影；

clear all

P = phantom('Modified Shepp-Logan',256)

R=radon(P)

figureimshow(R,[])

figure

imshow(P,[])title('仿体图')

%直接反投影法

l = pow2(nextpow2(size(R,1))-1)%重构图像的大小

P_1 = zeros(l,l)%用于存放重构后的图像

for i=1 : size(R,2)

    tmp = imrotate( repmat(R(:,i),1,size(R,1)),i-1,'bilinear' )

    tmp = tmp(floor(size(tmp,1)/2-l/2)+1:floor(size(tmp,1)/2+l/2),floor(size(tmp,2)/2-l/2)+1:floor(size(tmp,2)/2+l/2))

    P_1=P_1+tmp

end

P_1=P_1/size(R,2)

P_1=rot90(P_1)

figureimshow(P_1,[])title('直接反投影法')

%滤波反投影法

N=180

%滤波

H=size(R,1)

h=zeros((H*2-1),1)

for i=0:H-1

    if i==0

        h(H-i)=1/4

    elseif rem(i,2)==0

        h(H-i)=0

  姿御      h(H+i)=0

    else

        h(H-i)=-1/(i*pi)^2

        h(H+i)=-1/(i*pi)^2

  迹昌岩  end

end

x=zeros(H,N)

for i=1:N

    s=R(:,i)

    xx=conv(s',h')

    x(:,i)=xx(H:2*H-1)

end

%反投影

P_3=zeros(l,l)

for i=1:l

    for j=1:l

        for k=1:180

            theta=k/180*pi

            t=(j-l/2-0.5)*cos(theta)+(l/2+0.5-i)*sin(theta)+(H+1)/2

            t1=floor(t)

            t2=floor(t+1) 

            P_3(i,j)=P_3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k)

        end

    end

end

P_3=pi/N*P_3

figureimshow(P_3,[])title('滤波反投迅扮影法')

---