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('残差');