We consider this to be the most important part of the book. Each program was tested and has a comment statement to say which figures it generates; this is to put the reader at ease while running the program. Try to play with the programs by manipulating various parameters to see their effects. We have repeated the relevant chapter summary at the start of Section A.1 to A.6.
In this chapter we discussed the types of signals and the transformations a signal goes through before it is in a form suitable for input to a DSP processor. We introduced some common transforms, time and frequency domain interchanges, and the concepts of windowing, sampling and quantisation plus a model of an A/D converter. The last few sections give a flavour of DSP.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 1.1, 1.2 % Generations of Signals for 1st Chapter % 25th July 2003 clear;clf; % Narrow Band Signal fc=0.3; wn_upper=fc+fc*0.01; wn_lower=fc-fc*0.01; wn=[wn_lower wn_upper]; [b,a]=butter(2,wn); t=1:1000; u=randn(size(t)); y=filter(b,a,u); [sx,f]=spk(y,1); subplot(221);plot(t,y);grid; title(‘A’) xlabel(‘ Time ’) subplot(222);plot(f,sx);grid; title(‘B’) xlabel(‘ Frequency ’) %title (‘Fig 1.1’); pause; %print –depsc f1_1; %title(‘Fig 1.2’); %print –depsc f1_2; pause; % Broad band signal fc=0.6; wn_upper=fc+fc*0.1; wn_lower=fc-fc*0.1; wn=[wn_lower wn_upper]; [b,a]=butter(2,wn); t=1:1000; u=randn(size(t)); y=filter(b,a,u); [sx,f]=spk(y,1); subplot(221);plot(t,y);grid; title(‘A’); xlabel(‘ Time ’) subplot(222);plot(f,sx);grid; title(‘B’); xlabel(‘ Frequency ’) pause; %print –depsc f1_3; %title(‘Fig 1.3’) %title (‘Fig 1.4’) %print –depsc f1_4
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 1.3 % Generations of Signals for 1st Chapter % 25th July 2003 clear;clf; % Narrow Multi Band Signals fc=0.3;wn_upper=fc+fc*0.005;wn_lower=fc-fc*0.005; wn=[wn_lower wn_upper]; [b,a]=butter(2,wn); t=1:1000; u=randn(size(t)); y1=filter(b,a,u); fc=0.1; wn_upper=fc+fc*0.1; wn_lower=fc-fc*0.1; wn=[wn_lower wn_upper]; [b,a]=butter(2,wn); t=1:1000; u=randn(size(t)); y2=filter(b,a,u); fc=0.475; y3=0.1*sin(2*pi*fc*t); y=y1 + y2 + y3; [sx,f]=spk(y,1); subplot(221);plot(t,y);grid; xlabel(‘Time’);title(‘A’); subplot(222);plot(f,sx);grid; xlabel(‘Frequency’);title(‘B’); pause %print –depsc f1_5 %title (‘Fig 1.5’) %title (‘Fig 1.6’) %print –depsc f1_6
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 1.4 % Generations of Signals for 1st Chapter % 25th July 2003 clear;clf; % Statistical Signal t=1:5000; u=randn(size(t)); [x,y]=pdf_y(u,50); yhat=((exp(−(x.*x)*0.5))/(sqrt(2*pi))); scale=sum(yhat); yhat=yhat/scale; subplot(221);plot(t,u);grid; xlabel(‘Time’);title(‘A’); subplot(222);plot(x,y,‘o’,x,y,x,yhat,‘.–’);grid; xlabel(‘ Probability Density Function ’);title(‘B’); pause; print –depsc f1_7; %title (‘Fig 1.7’) %title (‘Fig 1.8’) %print –depsc f1_8;
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 1.5 and 1.6 clear;close; fs=10e3;T=1/fs; t = 1:1000; fc = 0.3; Q = 100;f1 = fc − (fc/Q);f2 = fc + (fc/Q); uk = randn(size(t)); wn = [f1 f2]; [b,a] = butter(2,wn); [h,w] = freqz(b,a,360); m = abs(h);F = w / (2*pi); xk = filter(b,a,uk); t=t*T; subplot(211);plot(t,xk);grid; % title(‘ Fig 1.9 ’); print –depsc f1_9; pause; [sxk,f] = spk(xk,1); subplot(211);plot(f*fs,sxk);grid; % title(‘ Fig 1.10 ’); xlabel(‘ Frequency in Hz. ’); print –depsc f11_0; pause subplot(211);plot(F*fs,m);grid; % title(‘ Fig 1.11 ’); xlabel(‘ Frequency in Hz. ’); print –depsc f11_1
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 1.7,1.8,1.9 %This program needs spk.m, expand.m clear;close; fs=10e3;T=1/fs; f=0.001; t=linspace(0,1000,1000); a=0.3 x=sin(2*pi*f*t)+a*sin(2*pi*f*t*3); impulse=ones([140]);t=t*T; timp_train=expand(impulse,26); imp_train=timp_train(1:1000); x_smple=x.*imp_train; x2=round((round(1000*x_smple))/200); x2=x2/5; subplot(221);plot(t,x,t,imp_train);grid; xlabel(‘ Seconds ’); %title(‘ Sampling of Continious signal ’); %title (‘ Fig 1.12 ‘) ; %subplot(221) ;plot(t,x_smple) ;grid; % xlabel(‘ Seconds ’); title(‘A‘); subplot(222) plot(t,x_smple,t,x,‘–.’);grid xlabel(‘ Seconds ’) ;title( ‘B’); % print –depsc f11_4; pause % print –depsc f11_2; pause; % print –depsc f11_3;pause; clf subplot(221);plot(t,x2,t,x,‘–.’);grid; title(‘ A ’); %title(‘ Fig 1.18 ’); subplot(222);plot(t,x_smple–x2);grid title(‘ B ’); xlabel(‘ Seconds ’); [sx_smple,f1]=spk(x_smple,1); [sx,f]=spk(x,1); % print –depsc f11_8 pause; clf %subplot(221) %plot(t,x) ;grid %title (‘ Continious Signal ’) %title (‘ Fig 1.15 ’) %subplot (223) %plot (f (1:10) *f s, sx(1:10));grid %title(‘ Fig 1.16 ’) %xlabel(‘ frequency Hz ’) subplot(211) plot(f1*fs,sx_smple*25);grid %title(‘ Fig 1.15 ’) xlabel(‘ frequency Hz ’) pause % print –depsc f11_5 %zoom;
This chapter was a refresher on some relevant engineering topics. It should help you with the rest of the book. We considered the autocorrelation function, representation of linear systems, and noise and its propagation in linear systems. We also discussed the need to know about systems reliability. Most problems in real life are inverse problems, so we introduced the Moore–Penrose pseudo-inverse. We concluded by providing an insight into binary number systems.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 2.3,2.4 clear;clf delta=0.1;f=100; omega=2*pi*f; b=omega*omega; a=[12*delta*omega omega*omega]; Ts=1/(10*f); sysc=tf(b,a) sysd=c2d(sysc,Ts,‘tustin’); sysd1=c2d(sysc,2*Ts,‘tustin’); [b1 a1]=tfdata(sysd,‘v’); [b2 a2]=tfdata(sysdl,‘v’); [y,tc]=impulse(sysc);tfinal=length(tc); td=1:tfinal; x=td*0; x(1)=1; y1=f ilter(b1,a1,x); y2=filter(b2,a2,x); delt=1;n1=length(x);delf = (1/delt)/n1; Sy = fft(y2); M = abs(Sy)*delt; n=fix(length(M)/2); [Mx i]=max(M(1:n)); mark=zeros(size(M)); mark(i)=M(i); f1=1:n; f1=f1*delf; Rtow=real(ifft(fftshift(M))); Rtow=fftshift(Rtow); T=1:length(Rtow); T=T−fix(length(Rtow)/2); T=T*delt; subplot(222);plot(f1,M(1:n),f1,mark(1:n));grid;title(‘B’) ylabel(‘ Power Spectrum ’); %title ( ‘ Fourier transformof R(tow) ’) %title ( ‘ Fig 2.5 ’) %print –depsc f2_5;pause; subplot(221);plot(T,Rtow);grid;title(‘A’) ylabel(‘ Auto Correlation ’); %title ( ‘ Auto Correlation function R(tow) ’); %title ( ‘ Fig 2.4 ’); print –depsc f2_4;pause td=td−1; %subplot(211) ;plot(tc,y) ;grid %subplot(212) ;plot(td,y1,td,y2);grid;pause; %subplot(211);plot(td(1:60),y2(1:60),td(1:60),y2(1:60),‘o’);grid subplot(211);stem(td(1:30), y2(1:30));grid ylabel(’ Impulse Response ’);xlabel(‘ Time ’) %title(‘Fig 2.3’); print –depsc f2_3
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 2.6 clear;close; delta=0.1; f=100; omega=2*pi*f; b1=omega*omega; a1=[12*delta*omega omega*omega]; sysc=tf(b1,a1); Ts1=1/(4*f); Ts2=1/(50*f); imax=10; Ts3=linspace(Ts1,Ts2,imax); for i=1:imax Ts=Ts3(i); sysd=c2d(sysc,Ts,‘tustin’); [b a]=tfdata(sysd,‘v’); q=abs(roots(a)); g1(i)=b(1); r1(i)=q(1); q=angle(roots(a)); theta1(i)=q(1); end imax=20;jmax=8; delt=linspace(0.05,0.9,imax); Ts1=1/(0.5*f) ;Ts2=1/(100*f); Ts3=linspace(Ts2,Ts1,jmax); Ts0=linspace(Ts3(1),Ts3(3),jmax); Ts0(jmax+1:2*jmax−3)=Ts3(4:jmax); jmax=length(Ts0); for j=1:jmax Ts=Ts0(j); k=1; for i=1:imax delta=delt(i); f=100; omega=2*pi*f; b1=omega*omega; a1=[12*delta*omega omega*omega]; sysc=tf(b1,a1); sysd=c2d(sysc,Ts,‘tustin’); [b a]=tfdata(sysd,‘v’); q1=abs(roots(a)); g2(i)=b(1); r2(k)=q1(1); q2=angle(roots(a)); theta2(k)=q2(1); k=k+1;r2(k)=q1(2);theta2(k)=q2(2);k=k+1; end %p=polyfit(delt,r2,1); %r2hat=p(1)*delt+p(2); %m(j)=p(1); c(j)=p(2); Tss(j)=Ts; A(j,:)=delt;B(j,:)=r2;C(j,:)=theta2; %subplot(211) ;plot(delt,r2, ‘o’, delt,r2hat) ;grid %subplot (212) ;plot (delt,theta2* (180/pi), ‘ o’) ;grid %polar (theta2,r2, ‘ o’);grid %theta(j)=mean(theta2*(180/pi)); %Ts end subplot(111);polar(C,B,‘o’);grid; %[xp,yp]=ginput(1);text(xp,yp,‘ Fig 2.7 ’); %[xp,yp]=ginput(1) ;text(xp,yp,‘<–––Constant Sampling Rate ’); %[xp,yp]=ginput(1) ;text(xp,yp,‘ Sensivity of System ’); pause; NN=(Tss/Ts1); %NN=1./NN; markx=linspace(0,180,180);marky=markx*0;marky(30)=1; marky(120)=1; subplot(221) plot(NN,B,markx/180,marky,‘b’);grid;xlabel(‘Sampling Frequency’); ylabel(‘Radial poleposition’);title(‘A’); subplot(222) plot(NN,C*180/pi);grid;xlabel(‘Sampling Frequency’); ylabel(‘Pole angle’) ,title(‘B’); %print –depsc f2_7
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 2.10 clear;close; n_first=100;n_middle=200; i=1:n_first; x=exp(−0.05*i); y=x; for i=1:n_middle x(n_first+i)=x(n_first); end for i=1:n_first x(n_first+n_middle+i)=y(n_first+1−i); end k=1:length(x); z=x*0;z=z+0.001; x=x+0.01; x=x*0.6; p1=n_first;p2=(p1+n_middle);p3=length(x); subplot(211); plot(k(1:p1),x(1:p1),‘o’,k(p1+1:p2),x(p1+1:p2), ‘−.’,k(p2+1:p3),x(p2+1:p3),‘*’,k,z) pause %print –depsc f2_7b
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 2.9 clear; kmax=5000;bins=100;order=5; k=1:kmax; g=9.81; mu_v=50; sig_v=0.2; mu_th=pi/3; sig_th=pi/30; temp=randn(size(k)); v=temp*sig_v;v=v+mu_v; temp=randn(size(k)); th=temp*sig_th; th=th+mu_th; d=v.*v.*sin(2*th)/g; [N,X]=hist(d,bins); pdf=N/kmax; [P,S]=polyfit(X,pdf,order); pdf_hat=polyval(P,X); subplot(221) plot(k,v);grid;title(‘ A ’);ylabel(‘ V in m/sec ’); subplot(222) plot(k,th*180/pi);grid;title(‘ B ’);ylabel(‘ ThetaDeg ’); subplot(212) plot(X,pdf,‘.’,X,pdf_hat);grid;title(‘ C ’);ylabel(‘Probability’);
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 2.14 %This program needs pdf_y.m clear;close; delta=0.1;f=100; omega=2*pi*f; b=omega*omega; a=[12*delta*omega omega*omega]; Ts=1/(10*f); sysc=tf(b,a) sysd=c2d(sysc,Ts,‘tustin’); [b1 a1]=tfdata(sysd,‘v’); plot(impulse(sysd));grid; a=0.1;b=0.8; for i=1:6000 u(i)=weibrnd(a,b); end z=filter(b1,a1,u); [x,y]=pdf_y(u,30); [x1,y1]=pdf_y(z,20); pause; subplot(211);plot(u);grid; %title(‘ Fig 2.8 ’); %print –depsc f2_8;pause; subplot(211);plot(x,y,‘*’,x,y,‘–’,x1,y1,‘.–’,x1,y1,‘o’);grid; %title(‘ Fig 2.9 ’); [xp,yp]=ginput(1);text(xp,yp,‘<––– Weibull distribution ’); [xp,yp]=ginput(1);text(xp,yp,‘<––– Normal distribution ’); %print –depsc f2_9
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 2.12,2.13 %This program needs pdf_y.m % Generations of Signals for 2nd Chapter % 25th July 2003 clear;clf; % Statistical Signal t=1:10000;u=randn(size(t)); [x,y]=pdf_y(u,50); sigma=1; xx=(x.*x)/(2*sigma); yhat=exp(−xx)/(sqrt(2*pi*sigma)); scale=sum(yhat);yhat=yhat/scale; [su,f1]=spk(u,1); subplot(222);plot(x,y,‘o’,x,y,x,yhat,‘.–’);grid; ylabel(‘ Statistics of uk ’); title(‘ B ’)’ %title(‘Fig 2.10’) %print –depsc f21_0;pause; subplot(221); plot(f1,su);grid; %title(‘Fig 2.11’) ylabel(‘White Noise ’); title(‘ A ’)’ %print –depsc f21_0;pause; %print –depsc f21_1 pause delta=0.1; %delta=0.01; f=100; omega=2*pi*f; b=omega*omega; a=[12*delta*omega omega*omega]; Ts=1/(10*f); sysc=tf(b,a); sysd=c2d(sysc,Ts,‘tustin’); [b1 a1]=tfdata(sysd,‘v’); y1=filter(b1,a1,u); %pause; [x,y]=pdf_y(y1,50); [sy1,f1]=spk(y1,1); gain=mean(sy1(1:150)); N=length(f1);m=abs(freqz(b1,a1,N)); sigma=1.5; %sigma=12; xx=(x.*x)/(2*sigma); yhat=exp(−xx)/(sqrt(2*pi*sigma)); scale=sum(yhat);yhat=yhat/scale; subplot(222);plot(x,y,‘o’,x,y,x,yhat,‘.–’);grid; ylabel(‘ PDF of yk ’);title(‘ B ’); %title(‘Fig 2.12’) subplot(221);plot(f1,sy1,‘.’,f1,m*gain,‘.–’);grid; ylabel(‘ Coloured Noise yk ’);title(‘ A ’); %title(‘Fig 2.13’) %xlabel(‘ Coloured Gaussian Noise – Ouput’) %print –depsc f21_2;pause; %print –depsc f21_3
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 2.17,2.18 clear;close sigma_x=0.2; sigma_y=0.2; x_val=−5;y_val=7; row_max=8; A=fix((rand([row_max 2])−0.5)*20); x1=x_val−randn( [1 row_max])*sigma_x; x2=y_val−randn([1 row_max])*sigma_y; c1=A(:,1);p1=c1.*x1’; c2=A(:,2);p2=c2.*x2’; p=fix(p1+p2); x3=inv(A‘*A)*A’*p; ek=A*x3−p; sigma=std(ek); x=x_val−3:0.1:x_val+3; for i=1:row_max q=A(i,:); y=−q(1)/q(2)*x+p(i)/q(2); B(i,:)=x;C(i,:)=y; end N=length(x); i=i+1; x_h=linspace(min(min(B)),max(max(B)), N); x_v=ones(size(x_h))*y_val; y_v=linspace(min(min(C)),max(max(C)), N); y_h=ones(size(y_v))*x_val; theta=linspace(0,2*pi,N); x=x3(1)−sigma*cos(theta); y=x3(2)−sigma*sin(theta); subplot(211); plot(B ‘,C’,‘.’,y_h,y_v,x_h,x_v,x,y,‘o’,x,y);grid; % print –depsc f21_6 subplot(212) stem(ek);grid % print –depsc f21_7
This chapter described ways to specify a given filter. We presented a variety of popular filters and their characteristics. We looked at IIR filters with real coefficients and truncated IIR filters for use in FIR designs. We also described a bank of filters with complex coefficients pointing to discrete Fourier transforms. We gave a simple presentation of adaptive filters and inverse problems having practical significance. We demonstrated an application of BFSK demodulation and adaptive phase shifting. And we looked at a practical inverse problem from target tracking. We added an important dimension by giving Kalman filter applications. Then we concluded this big topic with an application from array processing.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.1 clear;close; n=6;rdb=0.4;wn=0.4; [b,a]=cheby1(n,rdb,wn); %[b,a]=butter(n,wn); N=360; fn=linspace(0,0.5,N); sysd=tf(b,a,1); hk=impulse(sysd); m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi; nmax=length(fn);n1=round(wn*nmax*1.1); line_1=ones(size(fn(1:n1)))*0.92; line_2=ones(size(fn(1:n1)))*1.02; part_1=fn(1: n1); subplot(211);plot(fn,m,part_1,line_1,part_1,line_2);grid; print –depsc f3_1 %temp1=abs(roots(a)) ;temp=angle(roots(a)); % [theta,i]=sort(temp);r=temp1(i);pause; %subplot(111) ;polar(theta,r, ‘o’) ;grid; %subplot(212) ;plot(fn,phi,fn,phi_new) ;grid; %[xp, yp]=ginput (1); text(xp,yp,‘O<–––Non–linear Phase of IIR Filter ’) %[xp,yp]=ginput(1); text(xp,yp,‘O<–––Linear Phase of FIR filter ’) %print –depsc f3_8
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.2 clear;clg k=1:200;f=0.02;x=sin(2*pi*f*k)+0.2*rand(size(k)); subplot(211);stem(k,x);grid
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.9 clear;close; xk=1;xk1=0; for i=1:20 xk2=xk1;xk1=xk;xk=1.8*xk1−xk2; y(i)=xk;x(i)=i; end subplot(212); stem(x,y);grid text(xp,yp,‘1.8000 2.2400 2.2320 1.7776 0.9677−0.0358 −1.0321’); text(xp,yp,‘−1.8220 −2.2475 −2.2235 −1.7548 −0.9351 0.0715 1.0639’); %[xp,yp] =ginput (1); text(xp,yp,‘1.8435 2.2544 2.2144 1.7315 0.9023 −0.1073’); %[xp,yp]=ginput(1) ;text(xp,yp,‘ Fig 3.1 ’) print –depsc f3_2
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.10 clear;close; r=0.95;f=0.2; a=[1−2*r*cos(2*pi*f) r*r] ; b=[10−1];b=b*(1−r*r)*0.5;N=360; fn=linspace(0,0.5,N); sysd=tf(b,a,1); m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi;phi=phi/180; subplot(211);plot(fn,m,fn,phi);grid; [xp,yp]=ginput(1);text(xp,yp,‘<–––Phase +90 to −90 deg ’) [xp,yp]=ginput(1);text(xp,yp,‘<–––Amplitude Spectrrum ’) print –depsc f3_3
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 3.11,3.12 clear;close; N=50; fn=linspace(0,0.5,N); rmin=0.3;rmax=0.95; R=linspace(rmin,rmax,N); for k=1:N r=R(k); f=0.2; p=2*r*cos(2*pi*f); a1=[1−r*pr*r]; b1=[r*r −r*p 1]; phi=angle(freqz(b1,a1,N))*180/pi; i=find(phi > 0 );phi(i)=phi(i)− 360; A(k,:)=phi’; B(k,:)=fn; end subplot(211); plot(B,A,B(1,:),A(1,:),B(1,:),A(1,:),‘o’, B(N,:),A(N,:),B(N,:),A(N,:),‘d’);grid; [xp,yp]=ginput(1);text(xp,yp,‘<–––r = 0.3’); [xp,yp]=ginput(1);text(xp,yp,‘<–––r = 0.95’); print –depsc f3_5 pause; rmin=0.3;rmax=0.95; R=linspace(rmin,rmax,N); for k=1:N r=R(k);f=0.2; p=2*r*cos(2*pi*f); a2=[1−r]; b2=[−r 1]; phi=angle(freqz(b2,a2,N))*180/pi; A(k,:)=phi‘;B(k,:)=fn; end subplot(211) plot (B,A,B(1,:),A(1,:),B(1,:),A(1,:),‘o’, B(N,:),A(N,:),B(N,:),A(N,:),‘d’);grid; [xp,yp]=ginput(1);text(xp,yp,‘<–––r = 0.3’); [xp,yp]=ginput(1);text(xp,yp,‘<–––r = 0.95’); print –depsc f3_4
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.17 clear;close; r=0.95;f=0.2;p=2*cos(2*pi*f);g=(1−r*r)*0.5; a=[1 −r*p r*r]; b=[(1+r*r)/2−r*p (1+r*r)/2]; %b=b*(1−r*r)*0.5; N=360; fn=linspace(0,0.5,N); m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi; i=find( phi > 0 ); phi(i)=phi(i) − 360; phi1=((phi/180) + 1)*0.25+0.5; subplot(211);plot(fn,m,fn,phi1);grid; [xp,yp]=ginput(1);text(xp,yp,‘<–––Phase 0 to −360 deg ’) [xp,yp]=ginput(1);text(xp,yp,‘<–––Amplitude Spectrrum ’) print –depsc f3_6
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] clear;close; n=6;rdb=0.4;wn=0.4; [b,a]=cheby1(n,rdb,wn); %[b,a]=butter(n,wn); N=360; fn=linspace(0,0.5,N); sysd=tf(b,a,1);hk=impulse(sysd); m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi; nmax=length(fn);n1=round(wn*nmax*1.1); line_1=ones(size(fn(1:n1)))*0.92; line_2=ones(size(fn(1:n1)))*1.02; part_1=fn(1: n1); subplot(211);plot(fn,m,part_1,line_1,part_1,line_2);grid; temp1=abs(roots(a)); temp=angle(roots(a)); [theta,i]=sort(temp);r=temp1(i); subplot(111);polar(theta,r,‘o’);grid; x=r.*sin(theta);y=r.*cos(theta); [p,s]=polyfit(x,y,2);X=linspace(max(x), min(x),30); Y=polyval(p,X); theta_x=atan(X./Y);r_x=sqrt(X.*X+Y.*Y); hold; polar(theta_x,r_x,‘–.’); pause; %xlabel(‘Chebyshev Filter ’); %[xp,yp]=ginput (1) ;text (xp,yp,‘ Chebyshev Filter ’) print –depsc f31_4 %temp1=abs (roots (b)); temp=angle(roots(b)); [theta,i]=sort(temp);r=temp1(i); %pause; %subplot(111) ;polar(theta,r, ‘o’) ;grid; %subplot (212) ;plot (fn,phi,fn,phi_new) ;grid; %[xp,yp]=ginput(1) ;text(xp,yp, ‘O<–––Non–linear Phase of IIR Filter ’) %[xp,yp]=ginput(1) ;text(xp,yp,‘O<–––Linear Phase of FIR filter ’) %print –depsc f3_8
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.18,3.19 clear;close; n=6;rdb=0.4;wn=0.4; [b,a]=butter(n,wn); N=360; fn=linspace(0,0.5,N); sysd=tf(b,a,1);hk=impulse(sysd); m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi; nmax=length(fn);n1=round(wn*nmax*1.1); line_1=ones(size(fn(1:n1)))*0.92; line_2=ones(size(fn(1:n1)))*1.02; part_1=fn(1: n1); subplot(211);plot(fn,m,part_1,line_1,part_1,line_2);grid; temp1=abs(roots(a)); temp=angle(roots(a)); [theta,i]=sort(temp);r=temp1(i); subplot(111);polar(theta,r,‘o’);grid; x=r.*sin(theta);y=r.*cos(theta); [p,s]=polyfit(x,y,2);X=linspace(max(x), min(x),30); Y=polyval(p,X); theta_x=atan(X./Y);r_x=sqrt(X.*X+Y.*Y); hold; polar(theta_x,r_x,‘–.’); %xlabel(‘Butterworth Filter’); %[xp,yp]=ginput(1);text(xp,yp,‘Butterworth Filter ’) %[xp,yp]=ginput(1);text(xp,yp,‘<–––2nd Order pole trace ’) pause; %temp1=abs (roots (b)); temp=angle(roots(b));[theta,i]=sort(temp);r=temp1(i); %pause; %subplot(111) ;polar(theta,r, ‘o’) ;grid; %subplot (212) ;plot (fn,phi,fn,phi_new) ;grid; %[xp,yp]=ginput(1); %text(xp,yp, ‘O<–––Non–linear Phase of IIR Filter ’) %[xp, yp]=ginput(1); %text(xp,yp, ‘O<–––Linear Phase of FIR filter ’)
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.16 clear;close; t=1:500; f=0.15; x=sin(2*pi*f*t); phase=−30;phase=phase*pi/180; z=sin(2*pi*f*t+phase); plot(x,y);grid R=linspace(0.7,0.97,200); n=length(R); for i=1:n r=R(i); a=[1 −r];b=[1 −1/r]; y=filter(b,a,z); Jk(i)=sum((x−y).*(x−y)); %plot(t(1:20),x(1:20),‘o’,t(1:20),y(1:20),‘*’);grid; %pause; end subplot(211) plot(R,log(Jk)); xlabel(‘r’); ylabel(‘Jk(r)’);grid %pause %print –depsc f31_6a
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.3,3.4 clear;close; r=0.95;f=0.2; a=[1−2*r*cos(2*pi*f) r*r] ; b=[10−1];b=b*(1−r*r)*0.5;N=360; fn=linspace(0,0.5,N); sysd=tf(b,a,1);hk=impulse(sysd);nmax=8; m=abs(freqz(b,a,N)); phi=angle(freqz(b,a,N))*180/pi; bnew=hk(nmax:−1:1); bnew(nmax+1: 2*nmax−1)=hk(2:nmax); phi_new=angle(freqz(bnew,1,N))*180/pi; subplot(221); stem(bnew);grid;title(‘A’);ylabel(‘hk’); subplot(222);plot(fn,phi,fn,phi_new);grid;title(‘B’);ylabel(‘deg’); pause; %print –depsc f3_7; m_new=abs(freqz(bnew,1,N));m_new=m_new/(max(m_new)); subplot(211);plot(fn,m,fn,m_new);grid; pause; %print –depsc f3_8;pause;
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.21,3.22 clear close [m1,f]=freq(0.5); [m2,f]=freq(0.6); [m3,f]=freq(0.7); [m4,f]=freq(0.8); [m5,f]=freq(0.9); [m6,f]=freq(0.95); [m7,f]=freq(0.98); mn=[m1,m2,m3,m4,m5,m6,m7]; [m11,f]=freqt(36); [m21,f]=freqt(54); [m31,f]=freqt(72); [m41,f]=freqt(90); [m51,f]=freqt(108); [m61,f]=freqt(126); [m71,f]=freqt(144); m12=[m11 m21 m31 m41 m51 m61 m71]; subplot(211) plot(f,m12); xlabel (‘Normalised Frequency’); ylabel (‘Magnitude’) grid; %gtext(‘theta = 36’); %gtext(‘theta = 144’); pause %print –depsc f31_8 subplot(211) plot(f,mn); xlabel (‘Normalised Frequency’) ylabel(‘Magnitude’) grid; pause; %print –depsc f31_7 %gtext(‘r = 0.98’);
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 3.23,3.24 clear close range=2; N=200; delp=range /N; p=−1; r=0.95; r1=0.5*(1–r*r); k=1:150; f=0.2; power=0.5; u=sin(2*pi*f*k); noise=power*randn(size(u)); u=u+0*noise; for i=1:N b1=[r10−r1]; b2=[0 2*r]; a=[1 −2*r*p r*r]; p=p+delp; fhat=acos(p)/(2*pi); x=filter(b1,a,u); sk=filter(b2,a,x); ek=x − u; amp(i) =mean(ek.*ek); % F(i)=p; F(i)=fhat; slope(i)=2*mean(ek.*sk); end subplot(211) plot(F,amp); grid; xlabel p ylabel J(p) %title (‘ Objective function for r=0.95 at 3dB SNR ’); pause %print –depsc f31_9 subplot(211) plot(F,slope); grid; xlabel p ylabel ( ‘slope of [J(p)]’) pause %print –depsc f32_0 %zoom
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.29 clear;clf; b1=[10.6]; a1=[1 −1.7 1.53 −0.648]; [h,w]=freqz(b1,a1,180) ;m1=abs(h);f=w/(2*pi); subplot(221); plot(f,m1);grid;title(‘A’);xlabel(‘Frequency’); b2=0.999317; a2=[1 −2.16988 2.41286 −1.47436 0.365412]; [h,w]=freqz(b2,a2,180);m2=abs(h);f=w/(2*pi); subplot(222); plot(f,m2,f,m1,‘.–’);grid;title(‘B’);xlabel(‘Frequency’) %pause; %print –depsc f31_9a
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.31 clear;close; vw=15;N=2000; yw_max=N;yw=0:yw_max; delta_y=2*pi/yw_max; xw=500*sin(delta_y*yw); delta_t=0.25; k=0:N;t=k*delta_t; x_0=10000;y_0=8000; v_T=23;theta=120*(pi/180); x_T=x_0+v_T*cos(theta)*t; y_T=y_0+v_T*sin(theta)*t; noise=0.003*randn(size(t)); beta_m=atan2((x_T−xw),(y_T−yw))+noise; sum1=ones( [4 4]); sum1=sum1*0; sum2=[0 0 0 0]‘; for i=0:N u=[cos(beta_m(i+1)) −sin(beta_m(i+1)) t(i+1)*cos(beta_m(i+1)) −t(i+1)*sin(beta) yk=xw(i+1)*cos(beta_m(i+1))−yw(i+1)*sin(beta_m(i+1)); sum1=sum1 + u‘*u; sum2=sum2 + u‘*yk; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T=linspace(t(1),t(N+1),10); X=inv(sum1)*sum2; xhat_0=X(1); yhat_0=X(2); vhat=sqrt(X(3)*X(3)+X(4)*X(4)); xhat_T=xhat_0+vhat*cos(theta)*T; yhat_T=yhat_0+vhat*sin(theta)*T; subplot(221); plot(xw/1000,yw/1000,x_T/1000,y_T/1000,xhat_T/1000,yhat_T/ 1000,‘.–’);grid; title(‘A’); xlabel(‘ Target & Watcher ’); ylabel(‘ Km ’) subplot(222);plot(t,beta_m*180/pi);grid title(‘B’);xlabel(‘seconds’);ylabel(‘ beta−m in deg ’); pause;
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.32 % generate a linear reggression graph clear;close; t=linspace(0,2*pi,100); x=cos(t); y=x+0.1*randn(size(t)); subplot(211); plot(t,x,t,y,‘o’);grid; xlabel(‘ Time ’); label(‘ angle ’); pause % print –depsc f6_5
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.34 clear;clf; nmax=2000; t=1:nmax+9; f=0.05; xreal(1:1000)=2*sin(2*pi*f*t(1:1000)); xreal(1001:2000)=2*sin(2*pi*f*t(1001:2000)); theta=2*pi*f; % plot(u);grid;pause; % plot(xreal);pause; noise=randn(size(xreal)); y=xreal+0.5*noise; a=[2*cos(theta) −1; 10]; b=[0 0]‘; u=0; c=[10]; x_prev=[.01. 01]‘; p_prev=[10^6 0 ; 0 10^6]; N=0; s=100; for i=1:s N=N+1; pred_x = a*x_prev+b*u; pred_p = a*p_prev*a‘; err = c*pred_x − y(i); Y(i)=y(i) arr(i)=err; xhat = pred_x − ((pred_p*c‘)/(1+c*pred_p*c’))*err;phat = pred_p−(pred_p*c‘*c*pred_p)/(1+c*pred_p*c‘);Xfilter(i)=xhat(1); x_prev=xhat; p_prev=phat; end k=1:N; %plot(k,xreal(1:N),k,y(1:N),k,Xfilter);grid; subplot(211) plot(k,Y,‘o’,k,Xfilter);grid; pause print –depsc f6_9
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 3.26,3.25 %This program needs rbn_sig.m and expand.m files. clear; close; %N=5000; % generate BFSK data xx=sign(randn([1200])); fsk=rbn_sig(150,xx); ii=find(fsk==0); fsk=fsk*0.27;fsk(ii)=fsk(ii)+0.23;N=length(fsk); kk=1:N; z=sin(2*pi*fsk.*kk); [sz,ff]=spk(z,1);qq=800; subplot(221); plot(kk(1:qq), z(1: qq)*0.2+1.1,kk(1:qq),fsk(1:qq)*3);grid; title(‘A’);xlabel(‘ Time subplot(222);plot(ff,sz) ;grid; title(‘B’);xlabel(‘ Frequency ’); pause; % load fsk_data % K=1:N; % z=sin(2*pi*f*K); noise=0.2236*randn(size(z)); z=z+0.1*noise; lambda=0.9; mu =0.002; r=0.9; vk=0; delp =0; N=5000; range_p = 4; delP = range_p/N; xk1=0; xk=0; sk1=0; sk=0; uk1=0; uk=0; P = 0; for k=1:N uk2=uk1; uk1=uk; uk=z(k); xk2=xk1; xk1=xk; xk=r*p*xk1−r*r*xk2+(uk−uk2)*(1−r*r)*0.5; sk2=sk1; sk1=sk; sk=r*p*sk1−r*r*sk2+r*xk; ek=xk−uk; vk = vk*(lambda) + (1−lambda)*(ek*ek); grad = 2*ek*sk1; m(k)=mu; if (k > 1) delp = mu*grad/(vk+0.0001); end p = p − delp; ref(k)=acos(0.5*p)/(2*pi); end subplot(211) plot(kk(1:N),fsk(1:N),kk(1:N),ref,‘–.’) grid xlabel ‘Sample Number’ ylabel ‘p’ %title ‘ BFSK Demodulated Signal ’ pause; %print –depsc f61_1
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.35 clear;clf; t=linspace(0,2*pi,360); k=linspace(0,2*pi,60); z=sin(t);y=sin(k)+0.1*randn(size(k)); plot(t,z,k,y,‘o’); % Use this and Paint and PPT to generate f61_4
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 3.28 clear;clg; N=77; f=1/N; nmax=2000; K=1:nmax; phi=pi/20; delay=40; dly=1:delay; dly=dly*0; bufer=1: delay+1;bufer=bufer*0; Uk=cos(2*pi*f*K + phi ); rk=0.2; uk0=0; xk0=0; pk=0; noise=randn(size(Uk)); Uk=Uk+0.0*noise; loopGain=0.008;shift=0; for k=1:nmax; uk1=uk0; xk1=xk0; uk0=Uk(k); xk0=rk*xk1 − rk*uk0 + uk1; Xk(k)=xk0; zk=uk0*xk0; bufer(2:delay+1)=dly;bufer(1)=zk; dly=bufer(1:delay); ek=(dly(delay)−dly(1))/(delay); pk=pk+ek; Pk(k)=pk; rk=rk − loopGain*(pk+shift); Rk(k)=rk; T(k)=k; if abs(rk) > 0.95 rk=0.95; end end subplot(221) plot(Xk,Uk, ‘. ’);grid; xlabel(‘ xk ’); ylabel(‘ uk ’);title(‘ A ’); subplot(222) plot(T,Rk) ;grid;ylabel(‘ rk ’); title(‘ B ’);
In this chapter we used a matrix approach to explain the principles of the FFT. We did not used butterfly diagrams. Circular convolution was explained in depth. We presented a hardware scheme for implementation in real time. We looked at a problem of estimating frequency using a DFT. We elaborated a hardware structure for real-time implementation of continuous spectrum update. We ended by covering the principles of the network analyser, an example from RF systems.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 4.10 %This program needs p2f1.m, pulse.m and forier.m % demo package written by Dr K.V.Rangarao % this computes theoritical % Fourier Series Coefficients of a rectangular % pulse and compares with the FFT generated % Fourier Series Coefficients. clear;clf tow=0.1;T=1.2; f1=1/T;A=1.3;mm=9 N=2^mm; x=pulse(tow,T,N); x=x*A;delt=1e−3; n1=round(N/2); [f,a,b]=forier(x,delt); [a1,b1,F]=p2f 1(tow,f 1,A,N); %subplot (211) ;plot(x) ;grid; [aa ii]=max(a); s=ii−20;e=ii+20; subplot(221);plot(f(s:e),a(s:e),‘o’,f(s:e),a1(s:e),‘.–’); grid;title(‘A’) ylabel(‘ak’);xlabel(‘Co–sine’); subplot(222);plot(f(s:e),b(s:e),‘o’,f(s:e),b1(s:e),‘.–’); grid;title(‘B’) ylabel(‘bk’);xlabel(‘sine’); %print –depsc f4_4
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 4.11,4,12 clear;close; nx=50;ny=15; kx=1:nx; ky=1:ny; f=0.05; x=sin(2*pi*f*kx); y=cos(2*pi*f*ky); z=conv(x,y);t=1:length(z); %subplot (211); %plot (t,z); grid m=nx−ny; pad=ones([1 m])*0; y(ny+1:nx)=pad; zhat=real(ifft(fft(x).*fft(y))); t1=1:length(zhat); subplot(221) plot(t,z,‘.–’,t1,zhat,‘o’);grid;nn=length(zhat); e=z; e(1: nn)=e(1:nn)−zhat; subplot(222) plot(t,e,‘.–’);grid pause; %print –depsc f4_5 N=nx+ny−1; m=N−ny; pad=ones([1m])*0; y(ny+1:N)=pad; m=N−nx; pad=ones([1 m])*0; x(nx+1:N)=pad; zhat=real(ifft(fft(x).*fft(y))); t1=1:length(zhat); subplot(211) plot(t,z,‘.–’,t1,zhat,‘o’);grid;nn=length(zhat) e=z; e(1: nn)=e(1:nn)−zhat; %subplot(222) %plot(t,e,‘.–’);grid %pause %print –depsc f4_6
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 4.18 clear clg subplot(111) f=0.01; a=[1 −0.89944 0.404496]; b=[0.126 0.252 0.126]; t=1:2^ 11; zi=[0 0]; %rand(‘normal’); Fmax=50; F=linspace(0.05,0.37,Fmax); Kmax=length(F); for k=1:Kmax; f=F(k); x=cos(2*pi*f*t); [y,zf]=filter(b,a,x,zi); zi=zf; [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 uk=cos(2*pi*f1(i)*k); % yk=mag_true(i)*cos(2*pi*f1(i)*k+phse_true(i)*pi/180); yk=mag(i)*cos(2*pi*f1(i)*k+phase(i)*pi/180); temp(4−k)=yk; temp(7−k)=uk; Uk=cos(2*pi*f1(i)*K); Yk=mag(i)*cos(2*pi*f1(i)*K+phase(i)*pi/180); end % subplot (211); stem (Uk); % subplot(212);stem(Yk); temp(2:3)=−temp(2:3); Y(i)=temp(1);A(i,:)=temp(2:6); % pause; end z=Y′; p(1:2)=a(2:3);p(3:5)=b; phat=inv(A′*A)*A′*z;
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 4.15,4.16 clear;clf; w=20; T=100; f=0.3;records=6; x=ones([1 T]); x(w+1:T)=0*x(w+1:T); for i=1:records s=(i−1)*T+1; e=s+T−1; z(s:e)=x; end t=1:length(z);noise=1.0*randn(size(t)); y=sin(2*pi*f*t)+noise; puls=z.*y; subplot(211);plot(t,puls);grid; %title(‘A’); s1=1;e1=s1+w−1; for i=1:records s=(i−1)*T+1; e=s+w−1; xnew(s1:e1)=puls(s:e); s1=e1+1;e1=s1+w−1; end [sxnew,fx]=spk(xnew,1); %subplot(222) ;plot(xnew) ;title(‘B’) ;grid; pause; %print –depsc f6_2 subplot(221);plot(fx,sxnew,‘.–’);grid;title(‘A’);xlabel(‘ Freq ’); [smax,i]=max(sxnew);s=4;e=7; y1=sxnew(i−s:i+e);f1=fx(i−s:i+e); f2=linspace(f1(1),f1(length(f1)),100); p=polyfit(f1,y1,3) ;y1hat=polyval(p,f2); subplot(222); stem(f1,y1);grid hold; %plot(f1,y1,‘.–’,f2,y1hat);grid;title(‘B’) plot(f2,y1hat,‘.’);grid;title(‘B’); grid;xlabel(‘ Freq ’); pause; %print –depsc f6_2a
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 4.3 clear;close;clf; range=12; F_u=linspace(0.01,0.22,range); for k=1:range tiks=150; f_smp=200; f_u=F_u(k); x=1:tiks; x=x*0; x(tiks)=1; for i=1:f_smp s=1+(i−1)*tiks; e=s+tiks−1; u(s:e)=x; end b=1;a=[1−1]; z=filter(b,a,u); nz=length(z); z(nz)=z(nz−1); z=(z/max(z))*0.5; %z=(z/max(z))*0.25; t=1:nz; % Frequency Sweep xk=cos(2*pi*z.*t); % Unknown signal zk=sin(2*pi*f_u*t)+1*randn( [1nz]); % Band pass filter f_c=0.25; r=0.98; p=2*cos(2*pi*f_c); g=(1−r*r)*0.5; a=[1−r*pr*r] ; b=[10−1]*g; %................ uk=xk.*zk; yk=filter(b,a,uk); for i=1:f_smp s=1+(i−1)*tiks; e=s+tiks−1; blk=yk(s:e); rsp(i)=mean(blk.*blk); end %F=linspace(0,0.5,f_smp); F=linspace(0,1,f_smp); mid=round(f_smp/2); [rsp1,i1]=max(rsp(1:mid)); [rsp2,i2]=max(rsp(mid+1:f_smp)); %[rsp1,i1]=max(rsp); pk=rsp*0;pk(i1)=rsp1; pk(i2+mid)=rsp2; d(k)=F(i2+mid)−F(i1) ; %d(k)=F(i1); subplot(221) %plot(F(1:mid),rsp(1:mid),‘.–’,F(1:mid),pk(1:mid));grid;title(‘A’); plot(F,rsp,‘.–’,F,pk);grid;title(‘A’); pause end P=polyfit(d,F_u,1); Fhat=P(1)*d+P(2); subplot(222); plot(d,F_u,‘o’,d,Fhat);grid;title(‘ B ’); %pause %print –depsc f61_5b
In this chapter we described various hardware elements that are frequently used to convert the given analogue signal into digital form. We implemented a second-order BPF in SIMULINK and in fixed-point hardware then compared the results. We looked at a hardware realisation scheme that led us to consider pipelining and its associated problems. We gave an in-depth discussion on various methods of pipelining FIR and IIR filters. Then we looked at a complex algorithm for estimating frequency based on time series data and implemented on the DSP5630X using C. To implement this algorithm using FPGAs, we chose a design using FIR filters.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 5.8 clear; r2_o=0; r1_o=0; R2_o=0; R1_o=0; for i=1:1000 if i == 1 u=1; else u=0; end R1_i=(R2_o+u); R2_i=−0.9875*R1_o+1.608*R1_i; R1_o=R1_i; R2_o=R2_i; rl_i=fix(fix((r2_o+2 ^ *16*u))/256); r2_i=−252*r1_o+411*r1_i; r1_o=r1_i; r2_o=r2_i; x(i)=r1_o;X(i)=R1_o; end t=1:length(x); mark=ones(size(t))*255; m=abs(fft(x));M=abs(fft(256*X)); e=length(M) k=1:e/2;k=k/e; scale=1/max(M); m1=scale*m(1:e/2);M1=scale*M(1:e/2); [a p]=max(M1); subplot(211) plot(k(p−50:p+50),m1(p−50:p+50),‘–o’,k(p–50:p+50), M1(p−50:p+50));grid xlabel(‘ Normalised Frequency ’); ylabel(‘Response ’); pause; print –depsc f5_7a
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 5.11 clear; i=1; for phi=20:1:150 r=0.9; theta=phi*(pi/180); p=2*cos(theta); a1=−r*p; a2=r*r; P1=[1 a1 a2] ; P2=[1 −a1 a1*a1−a2]; P3=conv(P1,P2); a3=r^3*(2*p−p^3); a4=r^4*(p*p−1); P4=[1 0 0 a3 a4]; P4−P3; R=abs(roots(P2)); r1(i)=R(1);r2(i)=R(2);Theta(i)=phi; %pause i=i+1; end subplot(211) plot(Theta,r1,‘*’,Theta,r2,‘o’);grid %pause %print –depsc f51_1
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 5.13 close;clear kmax=5000; theta=20; f=theta/360; t=1:kmax; z=sin(2*pi*f*t)+0.2*randn(size(t)); a1=0.1; a2=0.1; r=0.98; v=0.0; Hk=[1e−9 0 0 0 1e−90 0 0 1e−9]; delp_k_old=[0 0 0]′; delp_k= [0 0 0]′; xk=0; xk1=0.1; sk=0; sk1=0; sk2=0; xk_corkt=0; zk=0; zk1=0; del_a1=0; del_a2=0; old_del_a1=0; old_del_a2=0; alpha1=0.95; delx_k=0; pk= [a1 a2 v]′; yk=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%% for k=1:1000; Dz=[1 −a1 −a2]; R=roots(Dz); rho(1,k)=abs(R(1)); theta(1,k)=angle(R(1)); rho(2,k)=abs(R(2)); theta(2,k)=angle(R(2)); alpha=0.98*alpha1; yk1=yk; yk=z(k); xk2=xk1; xk1=xk; xk=a1*xk1+a2*xk2; sk3=sk2; sk2=sk1; sk1=sk; sk=a1*sk1+a2*sk2+yk; ek=yk−xk; Jk=ek^2−v; Obj(k)=Jk; gk=[2*ek*sk12*ek*sk2 1]′; Hk_new=Hk/alpha; delH_num=Hk_new*gk*gk′*Hk_new; delH_den=1+gk′*Hk_new*gk; delH=delH_num/delH_den; Hk=Hk_new−delH; old_delp_k=delp_k; P(:,k)=pk; pk=pk+delp_k; a1=pk(1); a2=pk(2); v=pk(3); Jk_hat=Jk−gk′*old_delp_k; delp_k=Hk*gk*Jk_hat; del_a1=delp_k(1); del_a2=delp_k(2); delx_k=sk1*del_a1+sk2*del_a2; xk_corkt=xk+delx_k; xk=xk_corkt; alpha1=alpha1*0.9995+0.0005; x(k)=xk; y(k)=yk; T(k)=k; %%% A(k)=a1; B(k)=a2; C(k)=v; end subplot(221);plot(T,0.5*x,T,y);grid;zoom; subplot(222);plot(T,A,T,B,T,C);grid subplot(212);polar (theta,rho,‘o’);grid
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 5.14 clear; clf; load f513.dat; x1=f513(:,1); y1=f513(:,2); y2=f513(:,3); t=1:length(y1); s=20; e=90; subplot(221) plot(t,x1); grid; title(‘ A ’);xlabel(‘ Samples ’);ylabel(‘ Input ’); subplot(222) plot(t(s:e),y1(s:e),‘.’,t(s:e),y2(s:e));grid; title(‘ B ’);xlabel(‘ Samples ’);ylabel(‘ Output ’);
# include <stdio.h> # include <stdlib.h> # include <math.h> # include “f512matfun_f.c” # include “f512dsp.c” void main (void) { // test variables int k; int kmax = 2000; float theta=20.0,f,yk,xk,pi; floatpk[3] ={0.1,0.1,0}; FILE *fp_write; fp_write = fopen(“f512.dat”, “w”); f=theta/360.0; pi= (atan(1.0))*4.0; //randomise(); //printf(“ start ”); for (k=0;k < kmax; k++) { yk=sin(2*pi*f*k) + 0.1*randn(); xk = m_rem(yk,pk); // printf(“ xk = %.3e yk = %.3e ”,xk,yk); // getch(); fprintf (fp_write ,“%.2f %.2f %.2f %.2f %.2f ”, pk[0], pk[1], pk[2], xk,yk); } //printf(“finish ”); fclose(fp_write); } void mat_p( float *,int,int); void mat_mul( float *, float *, float *); void mat_mul_col( float *, float *, float *); void mat_mul_row( float *, float *, float *); void mat_trps( float *, float *); voidmat_eql( float *, float *); void ary_eql( float *, float *); voidmat_unt( float *); voidmat_scl( float *, float); void ary_scl( float *, float); void mat_mul_row_col( float *, float *, float *); voidmat_add( float *, float *, float *); void ary_add( float *, float *, float *); float mat_dot ( float *, float *); float randn(void); voidmat_mul( float *some_A, float *some_B, float *some_C) { int max_row_A = 3, max_col_A = 3; int max_col_B = 3; int col,row,k,index_A,index_B; float val_A,val_B,sum; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_B; col++) { sum = 0; for (k = 0; k < max_col_A; k++) { index_A = max_col_A*row + k; index_B = max_col_B*k + col; val_A = *(some_A + index_A); val_B = *(some_B + index_B); sum = sum + val_A*val_B; } *(some_C + max_col_A*row + col) = sum; } } } void mat_mul_row_col( float *some_A, float *some_B, float *some_C) { int max_col_C = 3, max_row_C = 3; int col,row; float val_A,val_B,sum; for ( row = 0; row < max_row_C; row++) { for ( col = 0; col < max_col_C; col++) { val_A = *(some_A + row); val_B = *(some_B + col); sum = val_A*val_B; *(some_C + max_col_C*row + col) = sum; } } } void mat_mul_col( float *temp_A, float *temp_B, float *temp_C) { int max_row_A = 3, max_col_A = 3; int row,k,index_A,index_B; float val_A,val_B,sum; for ( row = 0; row < max_row_A; row++) { sum = 0 ; for (k = 0;k < max_col_A; k++) { index_A = max_col_A*row + k; index_B = k; val_A = *(temp_A + index_A); val_B = *(temp_B + index_B); sum = sum + val_A*val_B; } *(temp_C + row) = sum; } } void mat_mul_row( float *temp_A, float *temp_B, float *temp_C) { int max_row_A = 3, max_col_A = 3; int col,k,index_A,index_B; float val_A,val_B,sum; for ( col = 0; col < max_row_A; col++) { sum = 0; for (k = 0;k < max_col_A; k++) { index_A = max_col_A*k + col; index_B = k; val_A = *(temp_A + index_A); val_B = *(temp_B + index_B); sum = sum + val_A*val_B; } *(temp_C + col) = sum; } } void mat_trps( float *some_A, float *some_B) { int max_row_A = 3, max_col_A = 3; int max_col_B = 3; int col,row,index_A,index_B; float val_A; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_B; col++) { index_A = max_col_A*row + col; index_B = max_col_B*col + row; val_A = *(some_A + index_A); *(some_B + index_B) = val_A; } } } void mat_eql( float *some_A, float *some_B) { int max_row_A = 3, max_col_A = 3; int max_col_B = 3; int col,row,index_A,index_B; float val_A; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_B; col++) { index_A = max_col_A*row + col; index_B = max_col_B*row + col; val_A = *(some_A + index_A); *(some_B + index_B) = val_A; } } } void ary_eql( float *some_A, float *some_B) { int max_row_A = 3; int row; float val_A; for (row= 0; row < max_row_A; row++) { val_A = *(some_A + row); *(some_B + row) = val_A; } } void ary_scl( float *some_A, float some_x) { int max_row_A = 3; int row; float val_A; for ( row = 0; row < max_row_A; row++) { val_A = *(some_A + row); *(some_A + row) = val_A*some_x; } } void ary_add( float *some_A, float *some_B, float *some_C) { int max_row_A = 3; int row; float val_A,val_B; for ( row = 0; row < max_row_A; row++) { val_A = *(some_A + row); val_B = *(some_B + row); *(some_C + row) = val_A + val_B; } } void mat_add( float *some_A, float *some_B, float *some_C) { int max_row_A = 3, max_col_A = 3; int max_col_B = 3; int col,row,index_A,index_B; float val_A,val_B; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_B; col++) { index_A = max_col_A*row + col; index_B = max_col_B*row + col; val_A = *(some_A + index_A); val_B = *(some_B + index_B); *(some_C + index_A) = val_A + val_B; } } } void mat_scl( float *some_A, float b) { int max_row_A = 3, max_col_A = 3; int col,row,index_A; float val_A; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_A; col++) { index_A = max_col_A*row + col; val_A = *(some_A + index_A); *(some_A + index_A) = val_A*b; } } } void mat_unt( float *some_A) { int max_row_A = 3, max_col_A = 3; int col,row,index_A; float val_A; for ( row = 0; row < max_row_A; row++) { for ( col = 0; col < max_col_A; col++) { index_A = max_col_A*row + col; val_A = 0; if ( col == row ) val_A =1; *(some_A + index_A) = val_A; } } } float mat_dot( float *temp_A, float *temp_B) { int max_row_A = 3; int col,k,index_A,index_B; float val_A,val_B,sum; sum = 0.0; for ( col = 0; col < max_row_A; col++) { index_A = col; index_B = col; val_A = *(temp_A + index_A); val_B = *(temp_B + index_B); sum = sum + val_A*val_B; } return(sum); } float randn(void) { float intmax,udf,sum; int rnd,k; intmax = 64000; sum = 0.0; for (k = 0; k < 12; k++) { rnd = rand(); udf = rnd; udf = (udf/intmax)−0.5; sum = sum + udf; } return(sum); } void mat_p( float *some_mat, int max_row, int max_col) { int col,row,index; float val; for (row=0; row < max_row; row++) { for (col=0; col < max_col; col++) { index = max_col*row + col; val= *(some_mat + index); printf (“ %2.3e ”,val); } printf(“ ”); } printf(“ ”); } float m_rem(float , float *); float m_rem(float yk, float *pk) { static float xk=0, xk1=0.1, xk2; static float sk=0, sk1=0, sk2=0, sk3, xk_corkt=0; static floatHk[3] [3] ={{1e-9,0,0},{0, 1e-9,0},{0,0,1e-9}}; static float Hk_new[3][3],delH_num[3][3],delH[3][3],delH_den; static float ek=0,Jk=0, Jk_hat=0; static float gk[3] = {0,0,0}; static float delp_k[3]={0, 0, 0},old_delp_k[3]={0, 0, 0}; static float alpha,alpha1=0.95, delx_k=0; static float temp1[3][3],temp2[3] [3]; static float scalar; alpha=0.98*alpha1; xk2=xk1; xk1=xk; xk=pk[0]*xk1+pk[1]*xk2; sk3=sk2; sk2=sk1; sk1=sk; sk=pk[0]*sk1+pk[1]*sk2+yk; ek=yk-xk; Jk=ek*ek-pk[2]; gk[0]=2*ek*sk1; gk[1]=2*ek*sk2; gk[2]= 1.0; mat_eql(Hk,Hk_new); scalar = 1.0/alpha; mat_scl(Hk_new,scalar); mat_mul_row_col(gk,gk,temp1); mat_mul(Hk_new,temp1,temp2); mat_mul(temp2,Hk_new,delH_num); mat_mul_col(Hk_new,gk,temp1); delH_den = 1 + mat_dot(gk,temp1); mat_eql(delH_num,delH); scalar = (-1.0/delH_den); mat_scl(delH,scalar); mat_add(Hk_new,delH,Hk); ary_eql(delp_k,old_delp_k); ary_add(pk,delp_k,pk); Jk_hat=Jk - mat_dot(gk,old_delp_k) mat_mul_col(Hk,gk,temp1); ary_scl(temp1,Jk_hat); ary_eql(temp1,delp_k); delx_k=sk1*delp_k[0] + sk2*delp_k[1]; xk_corkt=xk+delx_k; xk=xk_corkt; alpha1=alpha1*0.9995 + 0.0005; return(xk) ; } %Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] close;clear ! f51_2main load f512.dat; mrem=f512; p1 = mrem(:,1); p2 = mrem(:,2); p3 = mrem(:,3); p4 = mrem(:,4); p5 = mrem(:,5); k_c_code = length(p1); %%%%%%%%%%%%%%%%%%%%%%%%%%%% for k=1:k_c_code; Dz=[1-p1(k) -p2(k)]; R=roots(Dz); rho(1,k)=abs(R(1)); theta(1,k)=angle(R(1)); rho(2,k)=abs(R(2)); theta(2,k)=angle(R(2)); x(k)=p4(k); % Output Signal y(k)=p5(k); % Input Signal A(k)=p1(k); B(k)=p2(k); C(k)=p3(k); T(k)=k; %%% end s=800; e=850; subplot(221);plot(T(s:e),x(s:e),‘.’,T(s:e),x(s:e),T(s:e), y(s:e));grid; xlabel( ‘ Iterations ’ ) subplot(222);plot(T,A,T,B,T,C);grid xlabel( ‘ Iterations ’ ); ylabel( ‘ Parameters ’); % subplot(212);polar(theta,rho, ‘o’) ;grid;
In this chapter we gave a detailed treatment of a specific problem in direction finding, showing the need to understand systems engineering in order to implement DSP algorithms. It is very difficult to isolate DSP problems from other areas of engineering. Practical problems tend to cover many areas but centre on one area that demands the principal focus.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figure 6.4 clear;close; N=50;psi_deg=30; k=1:N;f=0.1; theta=30;A_x=sin(theta*pi/180);A_y=cos(theta*pi/180); psi_x=fix(rand(size(k))*psi_deg)*pi/180; psi_y=fix(rand(size(k))*psi_deg)*pi/180; x=A_x*cos(2*pi*f*k+psi_x); y=A_y*cos(2*pi*f*k+psi_y); subplot(221); plot(x,y,‘-’);grid;title(‘A’);xlabel(‘Analog’); %ylabel(‘-Y-’); u=ones(size(x)); A(:,1)=x′;A(:,2)=u′; p=inv(A′*A)*A‘*y’; X=linspace(min(x),max(x),100); Y=p(1)*X+p(2); subplot(222); plot(x,y,‘o’,X,Y);grid;title(‘B’);xlabel(‘Digital’); pause print -depsc f6_9b %ylabel(‘-Y-’);
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %Generates figures 6.1,6.2 clear; rot=2; for i=1:3 rot=rot−1; theta=0:180; theta=theta*pi/180; r=(sin(theta)). ^0.25; phi=theta + rot*pi/4; %polar (phi, r);grid;pause; X(:,i)=phi′; Y(:,i)=r′; end polar(X,Y);grid;pause; %print -depsc f6_3 A=ones([271 3]);A=A*0; A(1:181,1)=Y(1:181,1); A(46:226,2)=Y(1:181,2); A(91:271,3)=Y(1:181,3); for i=1:271 s=3*i-2;e=s+2; x(s:e)=A(i,:); end test=A(70,:); test=test+0.01*randn(size(test)); for i=1:271; ek=A(i,:)-test; Jk(i)=log(ek*ek′); end subplot(211);plot(Jk);grid;
This section gives subroutines needed by programs in the other sections.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function y=expand(x,r) % function y=expand(x,r) % introduces r zeros between % two succesive samples % e:matlab oolbox f j=1:length(x); k=j*r−r+1; n=length(k); y=zeros([1 n*r]); y(k)=x; clear k; clear j; return;
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function[m,f]=freq(r); t=40*pi/180; b=[10-1]; a=[1-2*r*cos(t) r*r]; [h,w]=freqz(b,a,180); f=w/(2*pi); m=abs(h); p=angle(h); return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function[m,f]=freqt(th); r=0.9; t=th*pi/180; b=[10-1]; a=[1-2*r*cos(t) r*r]; [h,w]=freqz(b,a,180); f=w/(2*pi); m=abs(h); p=angle(h); plot(f,m) grid return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %This program needs rbn_sig.m, expand.m function yy=rbn_sig(N,x) % generates a random binary telegraph signal % with +1 or 0 % where N = pulse width (samples) % and x is the binary sequence with +1 and −1 % ( bi-polar is a must) %x=sign(randn( [1 n])); % usage function yy=rbn_sig(N,x) y=expand(x,N); clear x; k=find( y > 0 ); j=diff(k); i=find (j == N); khat=k(i+1); k=find( y < 0 ); j=diff(k); i=find (j == N); khat_n=k(i+1); y(khat)=y(khat)*0; y(khat_n)=y(khat_n)*0; b=1; a=[1 -1]; subplot(211) plot(y);pause yy=abs(filter(b,a,y)); subplot(212) plot(yy) clear y; clear k; clear khat_n; clear khat; clear j; clear t; clear a; clear b; clear i; return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function [a,b,F]=p2f1(tow,f,A,n) % function [a,b,F]=p2f (tow,f, A,n) % converts a given pulse into % fourier coefficients % ‘a’ is a vector of cosine coefficeints % ‘b’ is a vector of sine coefficents. % ‘F’ is a vector giving the corresponding freq. % ‘tow’ is the pulse width. % ‘f’ is the frequency or the recipiocal of time-period. % ‘A’ is the amplitude of the pulse. % ‘n’ desired number of coefficients. N=fix(n/2); t=-N:1:N-1; tc=pi*f*tow; F=t*f; ta=A/pi; a=ta*(sin(2*tc*t))./t; b=2*ta*((sin(tc*t)). ^2)./t; %plot(b);pause; i=find(isnan(a)); v1=a(i+1); v2=a(i-1); a(i)=(v1+v2)/2; i=find(isnan(b)); v1=b(i+1);v2=b(i-1); b(i)=(v1+v2)/2; clear v1; clear v2; clear i; return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] %This program needs pulse.m function [y,delt]=pulse(tow,T,n); % function [y,delt]=pulse(tow,T,n); % generates a pulse of pulse width % ‘tow’ sec and time period of ‘T’ sec % of ‘n’ samples in the array ‘y’ delt=T/n; n1=round(tow/delt); t1=1:n1; y=1:n; y=0*y; y(1:n1)=ones(size(t1)); return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function [f,a,b]=forier(x,delt); % usage function [f,a,b]=forier(x,delt); % This finds the fourier coefficents of the % given periodic signal ‘x‘ % in dir d:matlab oolbox fforier.m Nx=length(x); frac=log(Nx)/log(2); if ( (frac - fix(frac)) > 0) m=fix(frac)+1; Nx=2 ^m; end; A=fftshift(fft(x,Nx)); n=round(length(A)/2); a=real(A)/n; b=-imag(A)/n; % A=fft(x); n=round(length(A)/2); % a=(real(A(1:n)))/n; b=-(imag(A(1:n)))/n; f=linspace(−0.5,0.5,length(A)); f=f/delt; clear A; return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function [sx,f]=spk(x,delt) % usage function [sx,f]=spk(x,delt); % This finds the spectrum of the given % vector x % returns ‘sx’ spectrum and ‘f’ the frequency% in dir e:matlab oolbox fspk.m n=length(x); frac=log(n)/log(2); if ( (frac-fix(frac)) > 0) m=fix(frac)+1; N=2 ^m; else N=n; end; % w=hamming(n); % x1=w.*x′; sxtemp=fft(x,N); f=linspace(0,0.5,N/2); f=f/delt; sx=(abs(sxtemp(1:N/2)))/(n/2); return
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] function [X,Y]=pdf_y(y,n) % function [X,Y]=pdf_y(y,n) % Finds the pdf of y specifed % by n bins and returns the % range and pdf in X and Y % resectively % e:matlab oolbox f l=length(y); [N,X]=hist(y,n); Y=N/l; clear l; return;
These programs demonstarte some of the key concepts for understanding the mathematics.
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] clear;close; x=−10:0.2:10; t=1:length(x); min_x=8; ref=ones(size(x))*min_x; y=2*(x-min_x).*(x-min_x)+6; xk=1; yk=2*(xk-min_x)*(xk-min_x)+6; gk=4*(xk-min_x); step = 1; for i=1:length(x)+20 t(i)=i;ref(i)=min_x; xl(i)=0;xl_prv(i)=0;xr(i)=0;xr_prv(i)=0; xk1=xk;yk1=yk;gk1=gk; xk=xk - step*yk/gk; yk=2*(xk-min_x)*(xk-min_x)+6+0.2*randn; gk=4*(xk-min_x)+0.1*randn; x1(i)=xk; y1(i)=yk; change_in_g=sign(gk1)-sign(gk); if change_in_g == 2 xk=xk1; step=step/1.5; end if change_in_g == −2 xk=xk1; step=step/1.5; end g(i)=change_in_g; %if g(i) == 2 %xl(i)=xk; xl_prv(i)=xk1; %end %if g(i) == −2 % xr(i)=xk; xr_prv(i)=xk1; %end end subplot(211);stem(t,g);grid; subplot(212);plot(t,x1,t,ref);grid; %pause; %subplot(211);plot(t,xl_prv, ‘o’, t,xl, ‘*’, t,ref);grid %subplot(212);plot(t,xr_prv, ‘o’, t,xr, ‘*’, t,ref);grid
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] clear;close A=rand([3 3])−0.5; A=fix(A*20); B=A+A′; c=rand([3 1])−0.5; c=fix(c*20); delB=c*c'; Z=inv(B+delB);invB=inv(B); Zhat=invB - invB*c*c′*invB/(1+c′*invB*c); Z-Zhat
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] clear;clf ;k=1; for i=30:1:50 theta=i*(pi/180);r=0.6; p=2*cos(theta); a3=r ^ 3*(2*p-p ^ 3); a4=r ^ 4*(p*p-1); b2=r*p;b3=r*r*(p*p-1); a1=-r*p; a2=r*r; a11=[1 0 0 a3 a4]; a11=1; a22=1; b22=[1 b2 b3]; a33=[1 a1 a2]; b33=1; [h1 w]=freqz(b33,a33,180); w=w/(2*pi); [h2 w]=freqz(b11,a11,180); w=w/(2*pi); m1=abs (h1); m2=abs(h2); rr=roots(a11); R=abs(rr);Theta=angle(rr); A(:,k)=R;B(:,k)=Theta; %polar(Theta,R,‘o’);grid; %pause k=k+1; end polar(B,A,‘o’);grid; %plot(w,h1,‘.-’,w,h2);grid
%Digital Signal Processing:A Practitioner's Approach %Dr.Kaluri Venkata Ranga Rao, [email protected] close;clear; N=3; i=1:N; a=[sum(i.*i) sum(i) sum(i) N ]; b=inv(a); t=1:N; u = ones(size(t)); x= [4 5 6]; c=[t*x′ u*x′] ′; mc=b*c a1=1; Nmax = 20; b1 = 2*(1:Nmax); b2 = (Nmax+1)*(ones(size(b1))); t1=1:200; noise=randn(size(t1))*6; x1=sin(2*pi*0.01*t1); x1=5*t1 + 1.5 + noise; y1=filter(b1,a1,x1); y2=filter(b2,a1,x1); y=((y2-y1)*6)/((Nmax-1)*(Nmax+1)*Nmax); subplot(211); plot(t1,y,t1,x1);grid subplot(212); plot(t1,y);grid;
Digital Signal Processing: A Practitioner's Approach K. V. Rangarao and R. K. Mallik
© 2005 John Wiley & Sons, Ltd