压缩感知的大神们，帮帮忙啊，一维信号的MP算法重构，稍微改动了下矩阵，出错了

clc
clear all;
close all;
M = 64;   %观测信号长度
N = 256;  %稀疏信号长度
K = 7;   %稀疏度
f1=50;    %  信号频率1
f2=100;   %  信号频率2
f3=200;   %  信号频率3
f4=400;   %  信号频率4
fs=800;   %  采样频率
ts=1/fs;  %  采样间隔
Ts=1:N;   %  采样序列
x0=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)+0.1*cos(2*pi*f3*Ts*ts)+0.9*cos(2*pi*f4*Ts*ts);

%------  高斯感知矩阵Phi   -------------
D=ones(M,N);
D(:,1:2:N)=0;
Phi = sqrt(1/M) * D;

%Phi = sqrt(1/M) * randn(M,N);
for i = 1:N 
    Phi(:,i) = Phi(:,i) / norm(Phi(:,i));
end
%-------- 测量向量 y  ----------
y = Phi * x0.';

%% -----2. MP Reconstruction ------------
Psi=fft(eye(N,N))/sqrt(N);                        %  傅里叶正变换矩阵
T=Phi*Psi';    

%  恢复矩阵(测量矩阵*正交反变换矩阵)
x = zeros(N,1);       %x0的逼近信号x
times =  K;           %迭代次数 = 稀疏度
g = zeros(N,1);       %余量和感知矩阵内积
r  = y;               %余量初始化为y
tic
for n=0:times
    g = T' * r;
    [val,K] = max( abs(g) ) ;
    x(K,1) = x(K,1) + g(K,1);
    r = r - g(K,1) * T(:,K);
end
toc
x1=real(Psi'*x); 
disp('MP relative error=');
disp( norm(x1'-x0)/norm(x0) );
figure(1),
hold on;
plot(x1,'k.-');plot(x0,'r');  legend('MP重构信号','原始信号');
figure(2),  plot(real(r),'b');  legend('残差');