clear all
clf;
%四十道雷克子信号的叠加
d=50;
v1=1400;
t00=2;
fm=23;
A=1;
FS0=500;
for p=1:40;
tx1=sqrt(t00.*t00+((p*d).*(p*d))./(v1.*v1));
[t,x]=ricker2(FS0,fm,tx1,A);
X(p,:)=x;
hold on
figure(1);plot(t,X(p,:)+p);
hold on
w=awgn(X(p,:),0,'measured');
W(p,:)=w;
hold on;
figure(2);plot(t,W(p,:)+p);hold on;
%将含噪信号缩放
%double S=zeros(40,1001);
a=0.5;
b=0;
s=(a-b)*(W(p,:)-min(W(p,:)))/(max(W(p,:))-min(W(p,:)))+b;
%将缩放后的信号通过TFPF算法滤波
mu=1;
z1(1)=s(1);
t1=1:length(t);
for i=2:length(s)
z1(i)=z1(i-1)+s(i);
end
z=exp(j*2*pi*mu.*z1);%首先对信号进行FM调频
h=ones(9,1);
[z2,t0,f0]=tfrpwv(z',t1,length(z'),h,0);%然后对调频信号求Wigner-Ville分布
y1=max(z2)/mu;%估计解析信号的峰值
y2=z2/mu;
for i=1:length(t);
for k1=1:length(f0);
if y2(k1,i)==y1(i)
y3(i)=f0(length(z')+1-k1);
end
end
end
%Y3(i,:)=y3;
%将信号反缩放得到期望信号
y4=(y3-b)*(max(w)-min(w))/(a-b)+min(w);
Y4(p,:)=y4;
%figure(1);plot(t,x+i);hold on ;figure(2);plot(t,w+i);
hold on;
figure(3);
plot(t,Y4(p,:)+p);
end