N=256; %信号长度
f1=50; %信号频率1
f2=100;%信号频率2
fs=800;%采样频率
ts=1/fs;%采样间隔
Ts=1:N;%采样序列
x=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)
fftx=fft(x);
k=4;
n=length(x)
m=30;
%产生一个随机的稀疏信号
x=2*(rand(1,n)-0.5);
x_sparse=rand(1,n);
x(x_sparse<0.95)=0;
plot(ts,x);
A=randn(m,n);%高斯随机测量矩阵
y=A*fftx';%投影距离y
err=10^(-8);
s=zeros(size(A,2),1);%size（A，2）表示A的列数
r=y;L=[];Psi=[];%初始化余量、索引值集合、感知矩阵
tic
while(i<2*k)(norm(r,2)<err)%norm(r)求矩阵r的2范数
h=A'*r;%求相关系数
[h_new,h_index]=sort(abs(h),'descend');%sort():对abs（h）的每一列按照降序排列，其中h_new为排序之后的矩阵，h_index中保存排序后响应元素在原矩阵中的序号
index=h_index(1);%将相关系数的最大值对应的索引值存入集合index
L=[L*index]';%更新支撑集
Psi=A(:,L);
c=Psi\y;%最小二乘法进行信号逼近
r=y-Psi*c;%余量更新
i=i+1;%迭代
end
s(L)=c;
toc
figure(1)
hold on;
subplot(2,1,1)
plot(Ts,x)
grid on;
xlable('Time');
ylable('Amplitude');
figure(2)
hold on;
len_s=length(s)
subplot(2,1,2)
plot(Ts,s,'r');
grid on;
xlable('Time');
ylable('Amplitude');