MATLAB and C Programs - Digital Signal Processing: A Practitioner's Approach - page 158

[y,zf]=filter(b,a,x,zi);
l=length(x);
xx=randn(size(t));
y=y+0.5*xx;
win=hamming(l);
win=win';
sx=fft(x.*win);
sy=fft(y.*win);
mx=(abs(sx))/(l);
my=(abs(sy))/(l);
phsex=atan2(imag(sx),real(sx))*180/pi;
phsey=atan2(imag(sy),real(sy))*180/pi;
%phsex=(angle (sx)) *180/pi;
%phsey=(angle (sy)) *180/pi;
[mx1,i1]=max(mx);
[my1,i2]=max(my);
gain=my1/mx1;
%phx(k)=phsex(i1);
%phy(k)=phsey(i2);
phse3=phsey(i2)-phsex(i1);
mag(k)=gain;
phase(k)=phse3;
f1(k)=f;
%h1=samplz(b,a,f);
Z=exp(sqrt(-1)*2*pi*f);
h1=polyval(b,Z)./polyval(a,Z);
mag_true(k)=abs(h1);
phse_true(k)=(angle(h1))*(180/pi);
%phse_true(k)=(atan2(imag(h1),real(h1))) * (180/pi);
%f=f+0.04
end
%phz=phx-phy;
%save fdls.dat f1 mag mag_true phase phse_true
%plot (f1,phy,f1,phx,f1,phz); grid
%pause
subplot(221);
plot(f1,mag,'o',f1,mag_true,'.-');grid;title(' A ');
subplot(222);
plot(f1,phase,'o',f1,phse_true,'.-');grid;title(' B ');
% FDLS starts from here
imax=length(f1);
K=0:10;
for i=1:imax-8;
for k=3:-1:1