clc;
clear all;
close all;
a1=clock; %记录时间
t=0:127;
y=(1/sqrt(280))*exp(-pi*((t-200)/280).^2).*cos(18*t+300)+(1/sqrt(250))*exp(-pi*((t-50)/250).^2).*cos(12*t+120)+randn(1,128);
figure (1);
subplot(2,1,1)
plot(t,y);
iterative_number=70; %迭代次数
[a,N]=size(y);
a_base=2;
j_min=0;
j_max=log2(N);
u_base=1/2;
p_min=0;
v_base=pi;
k_min=0;
w_base=pi/6;
i_min=0;
i_max=12;
%匹配追踪
signal_reconstruct1=zeros(1,N);
signal_r1=y;
for n=1:iterative_number
[proj,scale,translation,freq,phase]=select_best(signal_r1,N,a_base,j_min,j_max,u_base,p_min,v_base,k_min,w_base,i_min,i_max);
t=0:N-1;
g=(1/sqrt(scale))*exp(-pi*((t-translation)/scale).^2).*cos(freq*t+phase);
g=g/sqrt(sum(g.*g));
signal_reconstruct1=signal_reconstruct1+proj*g;
signal_r1=signal_r1-proj*g;
%n,scale,translation,freq,phase %数值跟踪
Ex1(n)=sum(signal_r1.*signal_r1);
end
%正交匹配追踪
signal_reconstruct2=zeros(1,N);
signal_r2=y;
pp=zeros(1,N);
uu=zeros(iterative_number,N);
for n=1:iterative_number
[proj,scale,translation,freq,phase]=select_best(signal_r2,N,a_base,j_min,j_max,u_base,p_min,v_base,k_min,w_base,i_min,i_max);
t=0:N-1;
g=(1/sqrt(scale))*exp(-pi*((t-translation)/scale).^2).*cos(freq*t+phase);
gu=zeros(1,N);
if n==1
pp=g;
uu(n,1:N)=pp/sqrt(sum(pp.*pp));
else
for k=1:n-1
gu=gu+sum(g.*uu(k,1:N))*uu(k,1:N);
end
pp=g-gu;
uu(n,1:N)=pp/sqrt(sum(pp.*pp));
end
signal_proj2=sum(signal_r2.*uu(n,1:N));
signal_reconstruct2=signal_reconstruct2+signal_proj2*uu(n,1:N);
signal_r2=signal_r2-signal_proj2*uu(n,1:N);
Ex2(n)=sum(signal_r2.*signal_r2);
%n,scale,translation,freq,phase %数值跟踪
end
t=0:N-1;
figure (2);
subplot(2,1,1)
plot(t,signal_r1);
subplot(2,1,2)
plot(t,signal_r2);
figure (3);
subplot(2,1,1)
plot(t,signal_reconstruct1);
subplot(2,1,2)
plot(t,signal_reconstruct1);
figure (4);
n=1:iterative_number
plot(n,Ex1,n,Ex2),legend('MP','OMP');
time=etime(clock,a1);
