c5_2_1.m

dsp入门与实践一书的源代码

clf;							
clear;							
w(1)=0.95;
r(1)=0.80;
N=1024;							
for t=1:N						
  w(t+1)=0.99*w(t)+0.01;
  r(t+1)=0.99*r(t)+0.00995;
end

dt=1/300;						
tt=0.:dt:dt*N;						
dis=1*(rand(size(tt)));					
y=100.*sin(2*pi*51*tt)+8*dis;				
p0=100/mean(y.*y);					

m(1)=0;e(1)=y(1);n(1)=0;a(1)=0;ei(1)=y(1);		
ef(1)=ei(1);yf(1)=y(1);
p(1)=p0/(w(1)+n(1)*n(1)*p0);
m(2)=-y(1)+r(2)*ei(1);
e(2)=y(2)+y(1)-m(2)*a(1);
n(2)=-yf(1)+r(2)*ef(1);
p(2)=p(1)/(w(2)+n(2)*n(2)*p(1));
a(2)=a(1)+p(1)*n(1)*e(1);
ei(2)=y(2)+y(1)-r(2)*r(2)*ei(1)-m(2)*a(2);
ef(2)=ei(2)-r(2)*ef(1)*a(2);
yf(2)=y(2)-r(2)*yf(1)*a(2);

for t=3:N							
  m(t)=-y(t-1)+r(t)*ei(t-1);
  e(t)=y(t)+y(t-2)-r(t).*r(t)*ei(t-2)-m(t)*a(t-1);
  n(t)=-yf(t-1)+r(t)*ef(t-1);
  p(t)=p(t-1)/(w(t)+n(t).*n(t)*p(t-1));
  a(t)=a(t-1)+p(t)*n(t)*e(t);
  ei(t)=y(t)+y(t-2)-r(t).*r(t)*ei(t-2)-m(t)*a(t);
  ef(t)=ei(t)-r(t).*r(t)*ef(t-2)-r(t)*ef(t-1)*a(t);
  yf(t)=y(t)-r(t).*r(t)*yf(t-2)-r(t)*yf(t-1)*a(t);
end

freq=(acos(-a/2))/(2*pi*dt);					
disp(freq(90:N))						
plot(freq(10:N))