7  Signal Processing toolbox

7.1   %asn elliptic integral

[y]=%asn(x,m)

• x : upper limit of integral (x>0) (can be a vector)
• m : parameter of integral (0<m<1)
• y : value of the integral

 y = integral from 0 to x of [1/(((1-t*t)^(1/2))(1-m*t*t)^(1/2))]

m=0.8;z=%asn(1/sqrt(m),m);K=real(z);Ktilde=imag(z);
x2max=1/sqrt(m);
x1=0:0.05:1;x2=1:((x2max-1)/20):x2max;x3=x2max:0.05:10;
x=[x1,x2,x3];
y=%asn(x,m);
rect=[0,-Ktilde,1.1*K,2*Ktilde];
plot2d(real(y)',imag(y)',1,'011',' ',rect)
//
deff('y=f(t)','y=1/sqrt((1-t^2)*(1-m*t^2))');
intg(0,0.9,f)-%asn(0.9,m)  //Works for real case only!

7.2   %k Jacobi's complete elliptic integral

[K]=%k(m)

• m : parameter of the elliptic integral 0<m<1 (m can be a vector)
• K : value of the elliptic integral from 0 to 1 on the real axis

K = integral from 0 to 1 of [(1-t^2)(1-m*t^2)]**(-1/2)

m=0.4;
%asn(1,m)
%k(m)

7.3   %sn Jacobi 's elliptic function

[y]=%sn(x,m)

• x : a point inside the fundamental rectangle defined by the elliptic integral; x is a vector of complex numbers
• m : parameter of the elliptic integral (0<m<1)
• y : result

m=0.36;
K=%k(m);
P=4*K; //Real period
real_val=0:(P/50):P;
plot(real_val,real(%sn(real_val,m)))
xbasc();
KK=%k(1-m);
Ip=2*KK;
ima_val1=0:(Ip/50):KK-0.001;
ima_val2=(KK+0.05):(Ip/25):(Ip+KK);
z1=%sn(%i*ima_val1,m);z2=%sn(%i*ima_val2,m);
plot2d([ima_val1',ima_val2'],[imag(z1)',imag(z2)']);
xgrid(3)

7.4   Signal Signal manual description

• analpf : analog low-pass filter
• buttmag : squared magnitude response of a Butterworth filter
• casc : creates cascade realization of filter
• cheb1mag : square magnitude response of a type 1 Chebyshev filter
• cheb2mag : square magnitude response of a type 1 Chebyshev filter
• chepol : recursive implementation of Chebychev polynomial
• convol : convolution of 2 discrete series
• ell1 mag : squared magnitude of an elliptic filter
• eqfir : minimax multi-band, linear phase, FIR filter
• eqiir : design of iir filter
• faurre : optimal lqg filter.
• lindquis : optimal lqg filter lindquist algorithm
• ffilt : FIR low-pass,high-pass, band-pass, or stop-band filter
• filter : compute the filter model
• find_freq : parameter compatibility for elliptic filter design
• findm : for elliptic filter design
• frmag : magnitude of the frequency responses of FIR and IIR filters.
• fsfirlin : design of FIR, linear phase (frequency sampling technique)
• fwiir : optimum design of IIR filters in cascade realization,
• iir : designs an iir digital filter using analog filter designs.
• iirgroup : group delay of iir filter
• iirlp : Lp IIR filters optimization
• group : calculate the group delay of a digital filter
• optfir : optimal design of linear phase filters using linear programming
• remezb : minimax approximation of a frequency domain magnitude response.
• kalm : Kalman update and error variance
• lev : resolve the Yule-Walker equations :
• levin : solve recursively Toeplitz system (normal equations)
• srfaur : square-root algorithm for the algebraic Riccati equation.
• srkf : square-root Kalman filter algorithm
• sskf : steady-state Kalman filter
• system : generates the next observation given the old state
• trans : transformation of standardized low-pass filter into low-pass, high-pass, band-pass, stop-band.
• wfir : linear-phase windowed FIR low-pass, band-pass, high-pass, stop-band
• wiener : Wiener estimate (forward-backward Kalman filter formulation)
• wigner : time-frequency wigner spectrum of a signal.
• window : calculate symmetric window
• zpbutt : Butterworth analog filter
• zpch1 : poles of a type 1 Chebyshev analog filter
• zpch2 : poles and zeros of a type 2 Chebyshev analog filter
• zpell : poles and zeros of prototype lowpass elliptic filter

• corr : correlation coefficients
• cspect : spectral estimation using the modified periodogram method.
• czt : chirp z-transform algorithm
• intdec : change the sampling rate of a 1D or 2D signal
• mese : calculate the maximum entropy spectral estimate
• pspect : auto and cross-spectral estimate
• wigner : Wigner-Ville time/frequency spectral estimation

• dft : discrete Fourier transform
• fft : fast flourier transform
• hilb : Hilbert transform centred around the origin.
• hank : hankel matrix of the covariance sequence of a vector process
• mfft : fft for a multi-dimensional signal

• lattn,lattp : recursive solution of normal equations
• phc : State space realisation by the principal hankel component approximation method,
• rpem : identification by the recursive prediction error method

• lgfft : computes p = ceil (log_2(x))
• sinc : calculate the function sin(2*pi*fl*t)/(pi*t)
• sincd : calculates the function Sin(N*x)/Sin(x)
• %k : Jacobi's complete elliptic integral
• %asn : .TP the elliptic integral :
• %sn : Jacobi 's elliptic function with parameter m
• bilt : bilinear transform or biquadratic transform.
• jmat : permutes block rows or block columns of a matrix

7.5   analpf create analog low-pass filter

[hs,pols,zers,gain]=analpf(n,fdesign,rp,omega)

• n : positive integer : filter order
• fdesign : string : filter design method : 'butt' or 'cheb1' or 'cheb2' or 'ellip'
• rp : 2-vector of error values for cheb1, cheb2 and ellip filters where only rp(1) is used for cheb1 case, only rp(2) is used for cheb2 case, and rp(1) and rp(2) are both used for ellip case. 0<rp(1),rp(2)<1
  
  1
• - for cheb1 filters 1-rp(1)<ripple<1 in passband
• - for cheb2 filters 0<ripple<rp(2) in stopband
• - for ellip filters 1-rp(1)<ripple<1 in passband 0<ripple<rp(2) in stopband
  
  0
• omega : cut-off frequency of low-pass filter in Hertz
• hs : rational polynomial transfer function
• pols : poles of transfer function
• zers : zeros of transfer function
• gain : gain of transfer function

//Evaluate magnitude response of continuous-time system 
hs=analpf(4,'cheb1',[.1 0],5)
fr=0:.1:15;
hf=freq(hs(2),hs(3),%i*fr);
hm=abs(hf);
plot(fr,hm)

7.6   buttmag response of Butterworth filter

[h]=buttmag(order,omegac,sample)

• order : integer : filter order
• omegac : real : cut-off frequency in Hertz
• sample : vector of frequency where buttmag is evaluated
• h : Butterworth filter values at sample points

//squared magnitude response of Butterworth filter
h=buttmag(13,300,1:1000);
mag=20*log(h)'/log(10);
plot2d((1:1000)',mag,[2],"011"," ",[0,-180,1000,20])

7.7   casc cascade realization of filter from coefficients

[cels]=casc(x,z)

• x : (4xN)-matrix where each column is a cascade element, the first two column entries being the numerator coefficients and the second two column entries being the denominator coefficients
• z : string representing the cascade variable
• cels : resulting cascade representation

x=[1,2,3;4,5,6;7,8,9;10,11,12]
cels=casc(x,'z')

7.8   cepstrum cepstrum calculation

fresp = cepstrum(w,mag)

• w : positive real vector of frequencies (rad/sec)
• mag : real vector of magnitudes (same size as w)
• fresp : complex vector

w=0.1:0.1:5;mag=1+abs(sin(w));
fresp=cepstrum(w,mag);
plot2d([w',w'],[mag(:),abs(fresp)])

7.9   cheb1mag response of Chebyshev type 1 filter

[h2]=cheb1mag(n,omegac,epsilon,sample)

• n : integer : filter order
• omegac : real : cut-off frequency
• epsilon : real : ripple in pass band
• sample : vector of frequencies where cheb1mag is evaluated
• h2 : Chebyshev I filter values at sample points

//Chebyshev; ripple in the passband
n=13;epsilon=0.2;omegac=3;sample=0:0.05:10;
h=cheb1mag(n,omegac,epsilon,sample);
plot(sample,h,'frequencies','magnitude')

7.10   cheb2mag response of type 2 Chebyshev filter

[h2]=cheb2mag(n,omegar,A,sample)

• n : integer ; filter order
• omegar : real scalar : cut-off frequency
• A : attenuation in stop band
• sample : vector of frequencies where cheb2mag is evaluated
• h2 : vector of Chebyshev II filter values at sample points

//Chebyshev; ripple in the stopband
n=10;omegar=6;A=1/0.2;sample=0.0001:0.05:10;
h2=cheb2mag(n,omegar,A,sample);
plot(sample,log(h2)/log(10),'frequencies','magnitude in dB')
//Plotting of frequency edges
minval=(-maxi(-log(h2)))/log(10);
plot2d([omegar;omegar],[minval;0],[2],"000");
//Computation of the attenuation in dB at the stopband edge
attenuation=-log(A*A)/log(10);
plot2d(sample',attenuation*ones(sample)',[5],"000")

7.11   chepol Chebychev polynomial

[Tn]=chepol(n,var)

• n : integer : polynomial order
• var : string : polynomial variable
• Tn : polynomial in the variable var

chepol(4,'x')

7.12   convol convolution

[y]=convol(h,x)
[y,e1]=convol(h,x,e0)

• x,h :input sequences (h is a "short" sequence, x a "long" one)
• e0 : old tail to overlap add (not used in first call)
• y : output of convolution
• e1 : new tail to overlap add (not used in last call)

x=1:3;
h1=[1,0,0,0,0];h2=[0,1,0,0,0];h3=[0,0,1,0,0];
x1=convol(h1,x),x2=convol(h2,x),x3=convol(h3,x),
convol(h1+h2+h3,x)
p1=poly(x,'x','coeff')
p2=poly(h1+h2+h3,'x','coeff')
p1*p2

7.13   corr correlation, covariance

[cov,Mean]=corr(x,[y],nlags)
[cov,Mean]=corr('fft',xmacro,[ymacro],n,sect)

[w,xu]=corr('updt',x1,[y1],w0)
[w,xu]=corr('updt',x2,[y2],w,xu)
 ...
[wk]=corr('updt',xk,[yk],w,xu)

• x : a real vector
• y : a real vector, default value x.
• nlags : integer, number of correlation coefficients desired.
• xmacro : a scilab external (see below).
• ymacro : a scilab external (see below), default value xmacro
• n : an integer, total size of the sequence (see below).
• sect : size of sections of the sequence (see below).
• xi : a real vector
• yi : a real vector,default value xi.
• cov : real vector, the correlation coefficients
• Mean : real number or vector, the mean of x and if given y

                n - m 
                 ====
                 \\                                                 1
        cov(m) =  >        (x(k)  - xmean) (y(m+k)      - ymean) * ---
                 /                                                  n
                 ====
                 k = 1

• - a function of type [xx]=xmacro(sect,istart) which returns a vector xx of dimension nsect containing the part of the sequence with indices from istart to istart+sect-1.
• - a fortran subroutine which performs the same calculation. (See the source code of dgetx for an example). n = total size of the sequence. sect = size of sections of the sequence. sect must be a power of 2. cov has dimension sect. Calculation is performed by FFT.

[w,xu]=corr('updt',x1,[y1],w0)
[w,xu]=corr('updt',x2,[y2],w,xu)
 ...
wk=corr('updt',xk,[yk],w,xu)

w0 = 0*ones(1,2*nlags);
nlags = power of 2.

x=%pi/10:%pi/10:102.4*%pi;
rand('seed');rand('normal');
y=[.8*sin(x)+.8*sin(2*x)+rand(x);.8*sin(x)+.8*sin(1.99*x)+rand(x)];
c=[];
for j=1:2,for k=1:2,c=[c;corr(y(k,:),y(j,:),64)];end;end;
c=matrix(c,2,128);cov=[];
for j=1:64,cov=[cov;c(:,(j-1)*2+1:2*j)];end;
rand('unif')
//
rand('normal');x=rand(1,256);y=-x;
deff('[z]=xx(inc,is)','z=x(is:is+inc-1)');
deff('[z]=yy(inc,is)','z=y(is:is+inc-1)');
[c,mxy]=corr(x,y,32);
x=x-mxy(1)*ones(x);y=y-mxy(2)*ones(y);  //centring
c1=corr(x,y,32);c2=corr(x,32);
norm(c1+c2,1)
[c3,m3]=corr('fft',xx,yy,256,32);
norm(c1-c3,1)
[c4,m4]=corr('fft',xx,256,32);
norm(m3,1),norm(m4,1)
norm(c3-c1,1),norm(c4-c2,1)
x1=x(1:128);x2=x(129:256);
y1=y(1:128);y2=y(129:256);
w0=0*ones(1:64);   //32 coeffs
[w1,xu]=corr('u',x1,y1,w0);w2=corr('u',x2,y2,w1,xu);
zz=real(fft(w2,1))/256;c5=zz(1:32);
norm(c5-c1,1)
[w1,xu]=corr('u',x1,w0);w2=corr('u',x2,w1,xu);
zz=real(fft(w2,1))/256;c6=zz(1:32);
norm(c6-c2,1)
rand('unif')
// test for Fortran or C external 
//
deff('[y]=xmacro(sec,ist)','y=sin(ist:(ist+sec-1))');
x=xmacro(100,1);
[cc1,mm1]=corr(x,2^3);
[cc,mm]=corr('fft',xmacro,100,2^3);
[cc2,mm2]=corr('fft','corexx',100,2^3);
[maxi(abs(cc-cc1)),maxi(abs(mm-mm1)),maxi(abs(cc-cc2)),maxi(abs(mm-mm2))]

deff('[y]=ymacro(sec,ist)','y=cos(ist:(ist+sec-1))');
y=ymacro(100,1);
[cc1,mm1]=corr(x,y,2^3);
[cc,mm]=corr('fft',xmacro,ymacro,100,2^3);
[cc2,mm2]=corr('fft','corexx','corexy',100,2^3);
[maxi(abs(cc-cc1)),maxi(abs(mm-mm1)),maxi(abs(cc-cc2)),maxi(abs(mm-mm2))]

7.14   cspect spectral estimation (periodogram method)

[sm,cwp]=cspect(nlags,ntp,wtype,x,y,wpar)

• x : data if vector, amount of input data if scalar
• y : data if vector, amount of input data if scalar
• nlags : number of correlation lags (positive integer)
• ntp : number of transform points (positive integer)
• wtype : string : 're','tr','hm','hn','kr','ch' (window type)
• wpar : optional window parameters for wtype='kr', wpar>0 and for wtype='ch', 0 < wpar(1) < .5, wpar(2) > 0
• sm : power spectral estimate in the interval [0,1]
• cwp : calculated value of unspecified Chebyshev window parameter

rand('normal');rand('seed',0);
x=rand(1:1024-33+1);
//make low-pass filter with eqfir
nf=33;bedge=[0 .1;.125 .5];des=[1 0];wate=[1 1];
h=eqfir(nf,bedge,des,wate);
//filter white data to obtain colored data 
h1=[h 0*ones(1:maxi(size(x))-1)];
x1=[x 0*ones(1:maxi(size(h))-1)];
hf=fft(h1,-1);   xf=fft(x1,-1);yf=hf.*xf;y=real(fft(yf,1));
sm=cspect(100,200,'tr',y);
smsize=maxi(size(sm));fr=(1:smsize)/smsize;
plot(fr,log(sm))

7.15   czt chirp z-transform algorithm

[czx]=czt(x,m,w,phi,a,theta)

• x : input data sequence
• m : czt is evaluated at m points in z-plane
• w : magnitude multiplier
• phi : phase increment
• a : initial magnitude
• theta : initial phase
• czx : chirp z-transform output

a=.7*exp(%i*%pi/6);
[ffr,bds]=xgetech(); //preserve current context
rect=[-1.2,-1.2*sqrt(2),1.2,1.2*sqrt(2)];
t=2*%pi*(0:179)/179;xsetech([0,0,0.5,1]);
plot2d(sin(t)',cos(t)',[2],"012",' ',rect)
plot2d([0 real(a)]',[0 imag(a)]',[3],"000")
xsegs([-1.0,0;1.0,0],[0,-1.0;0,1.0])
w0=.93*exp(-%i*%pi/15);w=exp(-(0:9)*log(w0));z=a*w;
zr=real(z);zi=imag(z);
plot2d(zr',zi',[5],"000")
xsetech([0.5,0,0.5,1]);
plot2d(sin(t)',cos(t)',[2],"012",' ',rect)
plot2d([0 real(a)]',[0 imag(a)]',[-1],"000")
xsegs([-1.0,0;1.0,0],[0,-1.0;0,1.0])
w0=w0/(.93*.93);w=exp(-(0:9)*log(w0));z=a*w;
zr=real(z);zi=imag(z);
plot2d(zr',zi',[5],"000")
xsetech(ffr,bds); //restore context

7.16   dft discrete Fourier transform

[xf]=dft(x,flag);

• x : input vector
• flag : indicates dft (flag=-1) or idft (flag=1)
• xf : output vector

n=8;omega = exp(-2*%pi*%i/n);
j=0:n-1;F=omega.^(j'*j);  //Fourier matrix
x=1:8;x=x(:);
F*x
fft(x,-1)
dft(x,-1)
inv(F)*x
fft(x,1)
dft(x,1)

7.17   ell1mag magnitude of elliptic filter

[v]=ell1mag(eps,m1,z)

• eps : passband ripple=1/(1+eps^2)
• m1 : stopband ripple=1/(1+(eps^2)/m1)
• z : sample vector of values in the complex plane
• v : elliptic filter values at sample points

deff('[alpha,beta]=alpha_beta(n,m,m1)',...
'if 2*int(n/2)=n then, beta=K1; else, beta=0;end;...
alpha=%k(1-m1)/%k(1-m);')
epsilon=0.1;A=10;  //ripple parameters
m1=(epsilon*epsilon)/(A*A-1);n=5;omegac=6;
m=find_freq(epsilon,A,n);omegar = omegac/sqrt(m)
%k(1-m1)*%k(m)/(%k(m1)*%k(1-m))-n   //Check...
[alpha,beta]=alpha_beta(n,m,m1)
alpha*%asn(1,m)-n*%k(m1)      //Check
sample=0:0.01:20;
//Now we map the positive real axis into the contour...
z=alpha*%asn(sample/omegac,m)+beta*ones(sample);
plot(sample,ell1mag(epsilon,m1,z))

7.18   eqfir minimax approximation of FIR filter

[hn]=eqfir(nf,bedge,des,wate)

• nf : number of output filter points desired
• bedge : Mx2 matrix giving a pair of edges for each band
• des : M-vector giving desired magnitude for each band
• wate : M-vector giving relative weight of error in each band
• hn : output of linear-phase FIR filter coefficients

hn=eqfir(33,[0 .2;.25 .35;.4 .5],[0 1 0],[1 1 1]);
[hm,fr]=frmag(hn,256);
plot(fr,hm),

7.19   eqiir Design of iir filters

[cells,fact,zzeros,zpoles]=eqiir(ftype,approx,om,deltap,deltas)

• ftype : filter type ('lp','hp','sb','bp')
• approx : design approximation ('butt','cheb1','cheb2','ellip')
• om : 4-vector of cutoff frequencies (in radians) om=[om1,om2,om3,om4], 0 <= om1 <= om2 <= om3 <= om4 <= pi. When ftype='lp' or 'hp', om3 and om4 are not used and may be set to 0.
• deltap : ripple in the passband. 0<= deltap <=1
• deltas : ripple in the stopband. 0<= deltas <=1
• cells : realization of the filter as second order cells
• fact : normalization constant
• zzeros : zeros in the z-domain
• zpoles : poles in the z-domain

[cells,fact,zzeros,zpoles]=...
eqiir('lp','ellip',[2*%pi/10,4*%pi/10],0.02,0.001)
transfer=fact*poly(zzeros,'z')/poly(zpoles,'z')

7.20   faurre filter

[Pn,Rt,T]=faurre(n,H,F,G,r0)

• n : number of iterations.
• H, F, G : estimated triple from the covariance sequence of y.
• r0 : E(yk*yk')
• Pn : solution of the Riccati equation after n iterations.
• Rt, T : gain matrix of the filter.

7.21   ffilt coefficients of FIR low-pass

[x]=ffilt(ft,n,fl,fh)

• ft : filter type where ft can take the values
  
  1
• "lp" : for low-pass filter
• "hp" : for high-pass filter
• "bp" : for band-pass filter
• "sb" : for stop-band filter
  
  0
• n : integer (number of filter samples desired)
• fl : real (low frequency cut-off)
• fh : real (high frequency cut-off)
• x : vector of filter coefficients

7.22   fft fast Fourier transform.

[x]=fft(a,-1)
[x]=fft(a,1)
x=fft(a,-1,dim,incr)
x=fft(a,1,dim,incr)

• x : real or complex vector. Real or complex matrix (2-dim fft)
• a : real or complex vector.
• dim : integer
• incr : integer

 a1=fft(a ,-1,dim1,incr1)
 a2=fft(a1,-1,dim2,incr2) ...

a=[1;2;3];n=size(a,'*');
norm(1/n*exp(2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*a -fft(a,1))
norm(exp(-2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*a -fft(a,-1))  

7.23   filter modelling filter

[y,xt]=filter(n,F,H,Rt,T)

• n : number of computed points.
• F, H : relevant matrices of the Markovian model.
• Rt, T : gain matrices.
• y : output of the filter.
• xt : filter process.

7.24   find_freq parameter compatibility for elliptic filter design

[m]=find_freq(epsilon,A,n)

• epsilon : passband ripple
• A : stopband attenuation
• n : filter order
• m : frequency needed for construction of elliptic filter

7.25   findm for elliptic filter design

[m]=findm(chi)

7.26   frfit frequency response fit

sys=frfit(w,fresp,order)
[num,den]=frfit(w,fresp,order)
sys=frfit(w,fresp,order,weight)
[num,den]=frfit(w,fresp,order,weight)

• w : positive real vector of frequencies (Hz)
• fresp : complex vector of frequency responses (same size as w)
• order : integer (required order, degree of den)
• weight : positive real vector (default value ones(w)).
• num,den : stable polynomials

w=0.01:0.01:2;s=poly(0,'s');
G=syslin('c',2*(s^2+0.1*s+2), (s^2+s+1)*(s^2+0.3*s+1));
fresp=repfreq(G,w);
Gid=frfit(w,fresp,4);
frespfit=repfreq(Gid,w);
bode(w,[fresp;frespfit])

7.27   frmag magnitude of FIR and IIR filters

[xm,fr]=frmag(num[,den],npts)

• npts : integer (number of points in frequency response)
• xm : mvector of magnitude of frequency response at the points fr
• fr : points in the frequency domain where magnitude is evaluated
• num : if den is omitted vector coefficients/polynomial/rational polynomial of filter
• num : if den is given vector coefficients/polynomial of filter numerator
• den : vector coefficients/polynomial of filter denominator

7.28   fsfirlin design of FIR, linear phase filters, frequency sampling technique

[hst]=fsfirlin(hd,flag)

• hd : vector of desired frequency response samples
• flag : is equal to 1 or 2, according to the choice of type 1 or type 2 design
• hst : vector giving the approximated continuous response on a dense grid of frequencies

//
//Example of how to use the fsfirlin macro for the design 
//of an FIR filter by a frequency sampling technique.
//
//Two filters are designed : the first (response hst1) with 
//abrupt transitions from 0 to 1 between passbands and stop 
//bands; the second (response hst2) with one sample in each 
//transition band (amplitude 0.5) for smoothing.
//
hd=[zeros(1,15) ones(1,10) zeros(1,39)];//desired samples
hst1=fsfirlin(hd,1);//filter with no sample in the transition
hd(15)=.5;hd(26)=.5;//samples in the transition bands
hst2=fsfirlin(hd,1);//corresponding filter
pas=1/prod(size(hst1))*.5;
fg=0:pas:.5;//normalized frequencies grid
plot2d([1 1].*.fg(1:257)',[hst1' hst2']);
// 2nd example
hd=[0*ones(1,15) ones(1,10) 0*ones(1,39)];//desired samples
hst1=fsfirlin(hd,1);//filter with no sample in the transition
hd(15)=.5;hd(26)=.5;//samples in the transition bands
hst2=fsfirlin(hd,1);//corresponding filter
pas=1/prod(size(hst1))*.5;
fg=0:pas:.5;//normalized frequencies grid
n=prod(size(hst1))
plot(fg(1:n),hst1);
plot2d(fg(1:n)',hst2',[3],"000");

7.29   group group delay for digital filter

[tg,fr]=group(npts,a1i,a2i,b1i,b2i)

• npts : integer : number of points desired in calculation of group delay
• a1i : in coefficient, polynomial, rational polynomial, or cascade polynomial form this variable is the transfer function of the filter. In coefficient polynomial form this is a vector of coefficients (see below).
• a2i : in coeff poly form this is a vector of coeffs
• b1i : in coeff poly form this is a vector of coeffs
• b2i : in coeff poly form this is a vector of coeffs
• tg : values of group delay evaluated on the grid fr
• fr : grid of frequency values where group delay is evaluated

z=poly(0,'z');
h=z/(z-.5);
[tg,fr]=group(100,h);
plot(fr,tg)

7.30   hank covariance to hankel matrix

[hk]=hank(m,n,cov)

• m : number of bloc-rows
• n : number of bloc-columns
• cov : sequence of covariances; it must be given as :[R0 R1 R2...Rk]
• hk : computed hankel matrix

//Example of how to use the hank macro for 
//building a Hankel matrix from multidimensional 
//data (covariance or Markov parameters e.g.)
//
//This is used e.g. in the solution of normal equations
//by classical identification methods (Instrumental Variables e.g.)
//
//1)let's generate the multidimensional data under the form :
//  C=[c_0 c_1 c_2 .... c_n]
//where each bloc c_k is a d-dimensional matrix (e.g. the k-th correlation 
//of a d-dimensional stochastic process X(t) [c_k = E(X(t) X'(t+k)], ' 
//being the transposition in scilab)
//
//we take here d=2 and n=64
//
c=rand(2,2*64)
//
//generate the hankel matrix H (with 4 bloc-rows and 5 bloc-columns)
//from the data in c
//
H=hank(4,5,c);
//

7.31   hilb Hilbert transform

[xh]=hilb(n[,wtype][,par])

• n : odd integer : number of points in filter
• wtype : string : window type ('re','tr','hn','hm','kr','ch') (default ='re')
• par : window parameter for wtype='kr' or 'ch' default par=[0 0] see the function window for more help
• xh : Hilbert transform

plot(hilb(51))

7.32   iir iir digital filter

[hz]=iir(n,ftype,fdesign,frq,delta)

• n : filter order (pos. integer)
• ftype : string specifying the filter type 'lp','hp','bp','sb'
• fdesign : string specifying the analog filter design ='butt','cheb1','cheb2','ellip'
• frq : 2-vector of discrete cut-off frequencies (i.e., 0<frq<.5). For lp and hp filters only frq(1) is used. For bp and sb filters frq(1) is the lower cut-off frequency and frq(2) is the upper cut-off frequency
• delta : 2-vector of error values for cheb1, cheb2, and ellip filters where only delta(1) is used for cheb1 case, only delta(2) is used for cheb2 case, and delta(1) and delta(2) are both used for ellip case. 0<delta(1),delta(2)<1
  
  1
• - for cheb1 filters 1-delta(1)<ripple<1 in passband
• - for cheb2 filters 0<ripple<delta(2) in stopband
• - for ellip filters 1-delta(1)<ripple<1 in passband and 0<ripple<delta(2) in stopband
  
  0

hz=iir(3,'bp','ellip',[.15 .25],[.08 .03]);
[hzm,fr]=frmag(hz,256);
plot2d(fr',hzm')
xtitle('Discrete IIR filter band pass  0.15<fr<0.25 ',' ',' ');
q=poly(0,'q');     //to express the result in terms of the ...
hzd=horner(hz,1/q) //delay operator q=z^-1

7.33   iirgroup group delay Lp IIR filter optimization

[lt,grad]=iirgroup(p,r,theta,omega,wt,td)
[cout,grad,ind]=iirlp(x,ind,p,[flag],lambda,omega,ad,wa,td,wt)

• r : vector of the module of the poles and the zeros of the filters
• theta : vector of the argument of the poles and the zeros of the filters
• omega : frequencies where the filter specifications are given
• wt : weighting function for and the group delay
• td : desired group delay
• lt, grad : criterium and gradient values

7.34   iirlp Lp IIR filter optimization

[cost,grad,ind]=iirlp(x,ind,p,[flag],lambda,omega,ad,wa,td,wt)

• x : 1X2 vector of the module and argument of the poles and the zeros of the filters
• flag : string : 'a' for amplitude, 'gd' for group delay; default case for amplitude and group delay.
• omega : frequencies where the filter specifications are given
• wa,wt : weighting functions for the amplitude and the group delay
• lambda : weighting (with 1-lambda) of the costs ('a' and 'gd' for getting the global cost.
• ad, td : desired amplitude and group delay
• cost, grad : criterium and gradient values

7.35   intdec Changes sampling rate of a signal

[y]=intdec(x,lom)

• x : input sampled signal
• lom : For a 1D signal this is a scalar which gives the rate change. For a 2D signal this is a 2-Vector of sampling rate changes lom=(col rate change,row rate change)
• y : Output sampled signal

7.36   jmat row or column block permutation

[j]=jmat(n,m)

• n : number of block rows or block columns of the matrix
• m : size of the (square) blocks

7.37   kalm Kalman update

[x1,p1,x,p]=kalm(y,x0,p0,f,g,h,q,r)

• f,g,h : current system matrices
• q, r : covariance matrices of dynamics and observation noise
• x0,p0 : state estimate and error variance at t=0 based on data up to t=-1
• y : current observation Output from the function is:
• x1,p1 : updated estimate and error covariance at t=1 based on data up to t=0
• x : updated estimate and error covariance at t=0 based on data up to t=0

7.38   lattn recursive solution of normal equations

[la,lb]=lattn(n,p,cov)

• n : maximum order of the filter
• p : fixed dimension of the MA part. If p= -1, the algorithm reduces to the classical Levinson recursions.
• cov : matrix containing the Rk's (d*d matrices for a d-dimensional process).It must be given the following way

                         |  R0 |
                         |  R1 |
                     cov=|  R2 |
                         |  .  |
                         |  .  |
                         |Rnlag|
  

• la : list-type variable, giving the successively calculated polynomials (degree 1 to degree n),with coefficients Ak

                       |Rp+1 Rp+2 . . . . . Rp+n  |
                       |Rp   Rp+1 . . . . . Rp+n-1|
                       | .   Rp   . . . . .  .    |
                       |                          |
   |I -A1 -A2 . . .-An|| .    .   . . . . .  .    |=0
                       | .    .   . . . . .  .    |
                       | .    .   . . . . .  .    |
                       | .    .   . . . . . Rp+1  |
                       |Rp+1-n.   . . . . . Rp    |

7.39   lattp lattp

[la,lb]=lattp(n,p,cov)

7.40   lev Yule-Walker equations (Levinson's algorithm)

[ar,sigma2,rc]=lev(r)

• r : correlation coefficients
• ar : auto-Regressive model parameters
• sigma2 : scale constant
• rc : reflection coefficients

      |R(0)   R(1)   ... R(N-1)|| ar(1) | |sigma2|
      |R(1)   R(0)   ... R(N-2)|| ar(2) | |  0   |
      |  :      :    ...    :  ||   :   |=|  0   |
      |  :      :    ...    :  ||   :   | |  0   |
      |R(N-1) R(N-2) ...  R(0) ||ar(N-1)| |  0   |
where
       R(i)=r(i-1)

7.41   levin Toeplitz system solver by Levinson algorithm (multidimensional)

[la,sig]=levin(n,cov)

• n : maximum order of the filter
• cov : matrix containing the R_k (d x d matrices for a d-dimensional process). It must be given the following way :

                  |  R0  |
                  |  R1  |
                  |  R2  |
                  |  .   |
                  |  .   |
                  | Rnlag|
  

• la : list-type variable, giving the successively calculated Levinson polynomials (degree 1 to n), with coefficients Ak
• sig : list-type variable, giving the successive mean-square errors.

                |R1   R2   . . . Rn  |
                |R0   R1   . . . Rn-1|
                |R-1  R0   . . . Rn-2|
                | .    .   . . .  .  |
|I -A1 .. -An|  | .    .   . . .  .  | = 0
                | .    .   . . .  .  |
                | .    .   . . .  .  |
                |R2-n R3-n . . . R1  |
                |R1-n R2-n . . . R0  |

//We use the 'levin' macro for solving the normal equations 
//on two examples: a one-dimensional and a two-dimensional process.
//We need the covariance sequence of the stochastic process.
//This example may usefully be compared with the results from 
//the 'phc' macro (see the corresponding help and example in it)
//
//
//1) A one-dimensional process
//   -------------------------
//
//We generate the process defined by two sinusoids (1Hz and 2 Hz) 
//in additive Gaussian noise (this is the observed process); 
//the simulated process is sampled at 10 Hz (step 0.1 in t, underafter).
//
t1=0:.1:100;rand('normal');
y1=sin(2*%pi*t1)+sin(2*%pi*2*t1);y1=y1+rand(y1);plot(t1,y1);
//
//covariance of y1
//
nlag=128;
c1=corr(y1,nlag);
c1=c1';//c1 needs to be given columnwise (see the section PARAMETERS of this help)
//
//compute the filter for a maximum order of n=10
//la is a list-type variable each element of which 
//containing the filters of order ranging from 1 to n; (try varying n)
//in the d-dimensional case this is a matrix polynomial (square, d X d)
//sig gives, the same way, the mean-square error
//
n=15;
[la1,sig1]=levin(n,c1);
//
//verify that the roots of 'la' contain the 
//frequency spectrum of the observed process y
//(remember that y is sampled -in our example 
//at 10Hz (T=0.1s) so that we need to retrieve 
//the original frequencies (1Hz and 2 Hz) through 
//the log and correct scaling by the frequency sampling)
//we verify this for each filter order
//
for i=1:n, s1=roots(la1(i));s1=log(s1)/2/%pi/.1;
//
//now we get the estimated poles (sorted, positive ones only !)
//
s1=sort(imag(s1));s1=s1(1:i/2);end;
//
//the last two frequencies are the ones really present in the observed 
//process ---> the others are "artifacts" coming from the used model size.
//This is related to the rather difficult problem of order estimation.
//
//2) A 2-dimensional process 
//   -----------------------
//(4 frequencies 1, 2, 3, and 4 Hz, sampled at 0.1 Hz :
//   |y_1|        y_1=sin(2*Pi*t)+sin(2*Pi*2*t)+Gaussian noise
// y=|   | with : 
//   |y_2|        y_2=sin(2*Pi*3*t)+sin(2*Pi*4*t)+Gaussian noise
//
//
d=2;dt=0.1;
nlag=64;
t2=0:2*%pi*dt:100;
y2=[sin(t2)+sin(2*t2)+rand(t2);sin(3*t2)+sin(4*t2)+rand(t2)];
c2=[];
for j=1:2, for k=1:2, c2=[c2;corr(y2(k,:),y2(j,:),nlag)];end;end;
c2=matrix(c2,2,128);cov=[];
for j=1:64,cov=[cov;c2(:,(j-1)*d+1:j*d)];end;//covar. columnwise
c2=cov;
//
//in the multidimensional case, we have to compute the 
//roots of the determinant of the matrix polynomial 
//(easy in the 2-dimensional case but tricky if d>=3 !). 
//We just do that here for the maximum desired 
//filter order (n); mp is the matrix polynomial of degree n
//
[la2,sig2]=levin(n,c2);
mp=la2(n);determinant=mp(1,1)*mp(2,2)-mp(1,2)*mp(2,1);
s2=roots(determinant);s2=log(s2)/2/%pi/0.1;//same trick as above for 1D process
s2=sort(imag(s2));s2=s2(1:d*n/2);//just the positive ones !
//
//There the order estimation problem is seen to be much more difficult !
//many artifacts ! The 4 frequencies are in the estimated spectrum 
//but beneath many non relevant others.
//

7.42   lgfft utility for fft

[y]=lgfft(x)

• x : real or complex vector

7.43   lindquist Lindquist's algorithm

[Pn,Rt,T]=lindquist(n,H,F,G,r0)

• n : number of iterations.
• H, F, G : estimated triple from the covariance sequence of y.
• r0 : E(yk*yk')
• Pn : solution of the Riccati equation after n iterations.
• RtP Tt : gain matrices of the filter.

7.44   mese maximum entropy spectral estimation

[sm,fr]=mese(x [,npts]);

• x : Input sampled data sequence
• npts : Optional parameter giving number of points of fr and sm (default is 256)
• sm : Samples of spectral estimate on the frequency grid fr
• fr : npts equally spaced frequency samples in [0,.5)

7.45   mfft multi-dimensional fft

[xk]=mfft(x,flag,dim)

• x : x(i,j,k,...) input signal in the form of a row vector whose values are arranged so that the i index runs the quickest, followed by the j index, etc.
• flag : (-1) FFT or (1) inverse FFT
• dim : dimension vector which gives the number of values of x for each of its indices
• xk : output of multidimensional fft in same format as for x

     x=[x(1,1,1),x(2,1,1),x(3,1,1),
        x(1,2,1),x(2,2,1),x(3,2,1),
              x(1,1,2),x(2,1,2),x(3,1,2),
              x(1,2,2),x(2,2,2),x(3,2,2),
                    x(1,1,3),x(2,1,3),x(3,1,3),
                    x(1,2,3),x(2,2,3),x(3,2,3)]

7.46   mrfit frequency response fit

sys=mrfit(w,mag,order)
[num,den]=mrfit(w,mag,order)
sys=mrfit(w,mag,order,weight)
[num,den]=mrfit(w,mag,order,weight)

• w : positive real vector of frequencies (Hz)
• mag : real vector of frequency responses magnitude (same size as w)
• order : integer (required order, degree of den)
• weight : positive real vector (default value ones(w)).
• num,den : stable polynomials

w=0.01:0.01:2;s=poly(0,'s');
G=syslin('c',2*(s^2+0.1*s+2),(s^2+s+1)*(s^2+0.3*s+1)); // syslin('c',Num,Den);
fresp=repfreq(G,w);
mag=abs(fresp);
Gid=mrfit(w,mag,4);
frespfit=repfreq(Gid,w);
plot2d([w',w'],[mag(:),abs(frespfit(:))])

7.47   phc Markovian representation

[H,F,G]=phc(hk,d,r)

• hk : hankel matrix
• d : dimension of the observation
• r : desired dimension of the state vector for the approximated model
• H, F, G : relevant matrices of the Markovian model

//
//This example may usefully be compared with the results from 
//the 'levin' macro (see the corresponding help and example)
//
//We consider the process defined by two sinusoids (1Hz and 2 Hz) 
//in additive Gaussian noise (this is the observation); 
//the simulated process is sampled at 10 Hz.
//
t=0:.1:100;rand('normal');
y=sin(2*%pi*t)+sin(2*%pi*2*t);y=y+rand(y);plot(t,y)
//
//covariance of y
//
nlag=128;
c=corr(y,nlag);
//
//hankel matrix from the covariance sequence
//(we can choose to take more information from covariance
//by taking greater n and m; try it to compare the results !
//
n=20;m=20;
h=hank(n,m,c);
//
//compute the Markov representation (mh,mf,mg)
//We just take here a state dimension equal to 4 :
//this is the rather difficult problem of estimating the order !
//Try varying ns ! 
//(the observation dimension is here equal to one)
ns=4;
[mh,mf,mg]=phc(h,1,ns);
//
//verify that the spectrum of mf contains the 
//frequency spectrum of the observed process y
//(remember that y is sampled -in our example 
//at 10Hz (T=0.1s) so that we need 
//to retrieve the original frequencies through the log 
//and correct scaling by the frequency sampling)
//
s=spec(mf);s=log(s);
s=s/2/%pi/.1;
//
//now we get the estimated spectrum
imag(s),
//

7.48   pspect cross-spectral estimate between 2 series

[sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar)

• x : data if vector, amount of input data if scalar
• y : data if vector, amount of input data if scalar
• sec_step : offset of each data window
• sec_leng : length of each data window
• wtype : window type (re,tr,hm,hn,kr,ch)
• wpar : optional parameters for wtype='kr', wpar>0 for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
• sm : power spectral estimate in the interval [0,1]
• cwp : unspecified Chebyshev window parameter

rand('normal');rand('seed',0);
x=rand(1:1024-33+1);
//make low-pass filter with eqfir
nf=33;bedge=[0 .1;.125 .5];des=[1 0];wate=[1 1];
h=eqfir(nf,bedge,des,wate);
//filter white data to obtain colored data 
h1=[h 0*ones(1:maxi(size(x))-1)];
x1=[x 0*ones(1:maxi(size(h))-1)];
hf=fft(h1,-1);   xf=fft(x1,-1);yf=hf.*xf;y=real(fft(yf,1));
//plot magnitude of filter
//h2=[h 0*ones(1:968)];hf2=fft(h2,-1);hf2=real(hf2.*conj(hf2));
//hsize=maxi(size(hf2));fr=(1:hsize)/hsize;plot(fr,log(hf2));
//pspect example
sm=pspect(100,200,'tr',y);smsize=maxi(size(sm));fr=(1:smsize)/smsize;
plot(fr,log(sm));
rand('unif');

7.49   remez Remez's algorithm

[an]=remez(nc,fg,ds,wt)

• nc : integer, number of cosine functions
• fg,ds,wt : real vectors
• fg : grid of frequency points in [0,.5)
• ds : desired magnitude on grid fg
• wt : weighting function on error on grid fg

 h = sum[a(n)*cos(wn)]

                 hn(1nc-1)=an(nc-12)/2;
                 hn(nc)=an(1);
                 hn(nc+12*nc-1)=an(2nc)/2;

7.50   remezb Minimax approximation of magnitude response

[an]=remezb(nc,fg,ds,wt)

• nc : Number of cosine functions
• fg : Grid of frequency points in [0,.5)
• ds : Desired magnitude on grid fg
• wt : Weighting function on error on grid fg
• an : Cosine filter coefficients

         hn(1:nc-1)=an(nc:-1:2)/2;
         hn(nc)=an(1);
         hn(nc+1:2*nc-1)=an(2:nc)/2;

// Choose the number of cosine functions and create a dense grid 
// in [0,.24) and [.26,.5)
nc=21;ngrid=nc*16;
fg=.24*(0:ngrid/2-1)/(ngrid/2-1);
fg(ngrid/2+1:ngrid)=fg(1:ngrid/2)+.26*ones(1:ngrid/2);
// Specify a low pass filter magnitude for the desired response
ds(1:ngrid/2)=ones(1:ngrid/2);
ds(ngrid/2+1:ngrid)=zeros(1:ngrid/2);
// Specify a uniform weighting function
wt=ones(fg);
// Run remezb
an=remezb(nc,fg,ds,wt)
// Make a linear phase FIR filter 
hn(1:nc-1)=an(nc:-1:2)/2;
hn(nc)=an(1);
hn(nc+1:2*nc-1)=an(2:nc)/2;
// Plot the filter's magnitude response
plot(.5*(0:255)/256,frmag(hn,256));
//////////////
// Choose the number of cosine functions and create a dense grid in [0,.5)
nc=21; ngrid=nc*16;
fg=.5*(0:(ngrid-1))/ngrid;
// Specify a triangular shaped magnitude for the desired response
ds(1:ngrid/2)=(0:ngrid/2-1)/(ngrid/2-1);
ds(ngrid/2+1:ngrid)=ds(ngrid/2:-1:1);
// Specify a uniform weighting function
wt=ones(fg);
// Run remezb
an=remezb(nc,fg,ds,wt)
// Make a linear phase FIR filter 
hn(1:nc-1)=an(nc:-1:2)/2;
hn(nc)=an(1);
hn(nc+1:2*nc-1)=an(2:nc)/2;
// Plot the filter's magnitude response
plot(.5*(0:255)/256,frmag(hn,256));

7.51   rpem RPEM estimation

[w1,[v]]=rpem(w0,u0,y0,[lambda,[k,[c]]])

• a,b,c : a=[a(1),...,a(n)], b=[b(1),...,b(n)], c=[c(1),...,c(n)]
• w0 : list(theta,p,phi,psi,l) where:
  
  1
• theta : [a,b,c] is a real vector of order 3*n
• p : (3*n x 3*n) real matrix.
• phi,psi,l : real vector of dimension 3*n
  
  0 During the first call on can take:

  theta=phi=psi=l=0*ones(1,3*n). p=eye(3*n,3*n)
  

• u0 : real vector of inputs (arbitrary size) (if no input take u0=[ ]).
• y0 : vector of outputs (same dimension as u0 if u0 is not empty). (y0(1) is not used by rpem).

y(t)+a(1)*y(t-1)+...+a(n)*y(t-n)=
b(1)*u(t-1)+...+b(n)*u(t-n)+e(t)+c(1)*e(t-1)+...+c(n)*e(t-n)

• lambda : optional parameter (forgetting constant) choosed close to 1 as convergence occur:

lambda(t)=alfa*lambda(t-1)+beta 

k(t)=mu*k(t-1)+nu 

7.52   sinc samples of sinc function

[x]=sinc(n,fl)

• n : number of samples
• fl : cut-off frequency of the associated low-pass filter in Hertz.
• x : samples of the sinc function

plot(sinc(100,0.1))

7.53   sincd sinc function

[s]=sincd(n,flag)

• n : integer
• flag : if flag = 1 the function is centred around the origin; if flag = 2 the function is delayed by %pi/(2*n)
• s : vector of values of the function on a dense grid of frequencies

plot(sincd(10,1)) 

7.54   srfaur square-root algorithm

[p,s,t,l,rt,tt]=srfaur(h,f,g,r0,n,p,s,t,l)

• h, f, g : convenient matrices of the state-space model.
• r0 : E(yk*yk').
• n : number of iterations.
• p : estimate of the solution after n iterations.
• s, t, l : intermediate matrices for successive iterations;
• rt, tt : gain matrices of the filter model after n iterations.
• p, s, t, l : may be given as input if more than one recursion is desired (evaluation of intermediate values of p).

7.55   srkf square root Kalman filter

[x1,p1]=srkf(y,x0,p0,f,h,q,r)

• f, h : current system matrices
• q, r : covariance matrices of dynamics and observation noise
• x0, p0 : state estimate and error variance at t=0 based on data up to t=-1
• y : current observation Output from the function is
• x1, p1 : updated estimate and error covariance at t=1 based on data up to t=0

7.56   sskf steady-state Kalman filter

[xe,pe]=sskf(y,f,h,q,r,x0)

• y : data in form [y0,y1,...,yn], yk a column vector
• f : system matrix dim(NxN)
• h : observations matrix dim(MxN)
• q : dynamics noise matrix dim(NxN)
• r : observations noise matrix dim(MxM)
• x0 : initial state estimate
• xe : estimated state
• pe : steady-state error covariance

7.57   system observation update

[x1,y]=system(x0,f,g,h,q,r)

• x0 : input state vector
• f : system matrix
• g : input matrix
• h : Output matrix
• q : input noise covariance matrix
• r : output noise covariance matrix
• x1 : output state vector
• y : output observation

     x1=f*x0+g*u
     y=h*x0+v

7.58   trans low-pass to other filter transform

hzt=trans(pd,zd,gd,tr_type,frq)

• hz : input polynomial
• tr_type : type of transformation
• frq : frequency values
• hzt : output polynomial

7.59   wfir linear-phase FIR filters

[wft,wfm,fr]=wfir(ftype,forder,cfreq,wtype,fpar)

• ftype : string : 'lp','hp','bp','sb' (filter type)
• forder : Filter order (pos integer)(odd for ftype='hp' or 'sb')
• cfreq : 2-vector of cutoff frequencies (0<cfreq(1),cfreq(2)<.5) only cfreq(1) is used when ftype='lp' or 'hp'
• wtype : Window type ('re','tr','hm','hn','kr','ch')
• fpar : 2-vector of window parameters. Kaiser window fpar(1)>0 fpar(2)=0. Chebyshev window fpar(1)>0, fpar(2)<0 or fpar(1)<0, 0<fpar(2)<.5
• wft : time domain filter coefficients
• wfm : frequency domain filter response on the grid fr
• fr : Frequency grid

7.60   wiener Wiener estimate

[xs,ps,xf,pf]=wiener(y,x0,p0,f,g,h,q,r)

• f, g, h : system matrices in the interval [t0,tf]
  
  1
• f =[f0,f1,...,ff], and fk is a nxn matrix
• g =[g0,g1,...,gf], and gk is a nxn matrix
• h =[h0,h1,...,hf], and hk is a mxn matrix
  
  0
• q, r : covariance matrices of dynamics and observation noise
  
  1
• q =[q0,q1,...,qf], and qk is a nxn matrix
• r =[r0,r1,...,rf], and gk is a mxm matrix
  
  0
• x0, p0 : initial state estimate and error variance
• y : observations in the interval [t0,tf]. y=[y0,y1,...,yf], and yk is a column m-vector
• xs : Smoothed state estimate xs= [xs0,xs1,...,xsf], and xsk is a column n-vector
• ps : Error covariance of smoothed estimate ps=[p0,p1,...,pf], and pk is a nxn matrix
• xf : Filtered state estimate xf= [xf0,xf1,...,xff], and xfk is a column n-vector
• pf : Error covariance of filtered estimate pf=[p0,p1,...,pf], and pk is a nxn matrix

7.61   wigner 'time-frequency' wigner spectrum

[tab]=wigner(x,h,deltat,zp)

• tab : wigner spectrum (lines correspond to the time variable)
• x : analyzed signal
• h : data window
• deltat : analysis time increment (in samples)
• zp : length of FFT's. %pi/zp gives the frequency increment.

7.62   window symmetric window

[win_l,cwp]=window(wtype,n,par)

• wtype : window type (re,tr,hn,hm,kr,ch)
• n : window length
• par : parameter 2-vector (kaiser window: par(1)=beta>0) (Chebychev window par =[dp,df]), dp =main lobe width (0<dp<.5), df=side lobe height (df>0)
• win : window
• cwp : unspecified Chebyshev window parameter

7.63   yulewalk least-square filter design

Hz = yulewalk(N,frq,mag) 

• N : integer (order of desired filter)
• frq : real row vector (non-decreasing order), frequencies.
• mag : non negative real row vector (same size as frq), desired magnitudes.
• Hz : filter B(z)/A(z)

                  n-1         n-2            
      B(z)   b(1)z     + b(2)z    + .... + b(n)
H(z)= ---- = ---------------------------------
                n-1       n-2
      A(z)    z   + a(2)z    + .... + a(n)

f=[0,0.4,0.4,0.6,0.6,1];H=[0,0,1,1,0,0];Hz=yulewalk(8,f,H);
fs=1000;fhz = f*fs/2;  
xbasc(0);xset('window',0);plot2d(fhz',H');
xtitle('Desired Frequency Response (Magnitude)')
[frq,repf]=repfreq(Hz,0:0.001:0.5);
xbasc(1);xset('window',1);plot2d(fs*frq',abs(repf'));
xtitle('Obtained Frequency Response (Magnitude)')

7.64   zpbutt Butterworth analog filter

[pols,gain]=zpbutt(n,omegac)

• n : integer (filter order)
• omegac : real (cut-off frequency in Hertz)
• pols : resulting poles of filter
• gain : resulting gain of filter

7.65   zpch1 Chebyshev analog filter

[poles,gain]=zpch1(n,epsilon,omegac)

• n : integer (filter order)
• epsilon : real : ripple in the pass band (0<epsilon<1)
• omegac : real : cut-off frequency in Hertz
• poles : resulting filter poles
• gain : resulting filter gain

 H(s)=gain/poly(poles,'s')

7.66   zpch2 Chebyshev analog filter

[zeros,poles,gain]=zpch2(n,A,omegar)

• n : integer : filter order
• A : real : attenuation in stop band (A>1)
• omegar : real : cut-off frequency in Hertz
• zeros : resulting filter zeros
• poles : resulting filter poles
• gain : Resulting filter gain

H(s)=gain*poly(zeros,'s')/poly(poles,'s')

7.67   zpell lowpass elliptic filter

[zeros,poles,gain]=zpell(epsilon,A,omegac,omegar)

• epsilon : real : ripple of filter in pass band (0<epsilon<1)
• A : real : attenuation of filter in stop band (A>1)
• omegac : real : pass band cut-off frequency in Hertz
• omegar : real : stop band cut-off frequency in Hertz
• zeros : resulting zeros of filter
• poles : resulting poles of filter
• gain : resulting gain of filter

7.68   arma Scilab arma library

• armac : this function creates a description as a list of an ARMAX process A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t)
• armap : Display in the file out or on the screen the armax equation associated with ar
• armax : is used to identify the coefficients of a n-dimensional ARX process A(z^-1)y= B(z^-1)u + sig*e(t)
• armax1 : armax1 is used to identify the coefficients of a 1-dimensional ARX process A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t)
• arsimul : armax trajectory simulation
• arspec : Spectral power estimation of armax processes. Test of mese and arsimul
• exar1 : An Example of ARMAX identification ( K.J. Astrom) The armax process is described by : a=[1,-2.851,2.717,-0.865] b=[0,1,1,1] d=[1,0.7,0.2]
• exar2 : ARMAX example ( K.J. Astrom). A simulation of a bi dimensional version of the example of exar1.
• exar3 : Spectral power estimation of arma processes from Sawaragi et all where a value of m=18 is used. Test of mese and arsimul
• gbruit : noise generation
• narsimul : armax simulation ( using rtitr)
• odedi : Simple tests of ode and arsimul. Tests the option 'discret' of ode
• prbs_a : pseudo random binary sequences generation
• reglin : Linear regression

7.69   armac Scilab description of an armax process

[ar]=armac(a,b,d,ny,nu,sig)

• a=[Id,a1,..,a_r] : is a matrix of size (ny,r*ny)
• b=[b0,.....,b_s] : is a matrix of size (ny,(s+1)*nu)
• d=[Id,d1,..,d_p] : is a matrix of size (ny,p*ny);
• ny : dimension of the output y
• nu : dimension of the output u
• sig : a matrix of size (ny,ny)

   A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t)

 
a=[1,-2.851,2.717,-0.865].*.eye(2,2)
b=[0,1,1,1].*.[1;1];
d=[1,0.7,0.2].*.eye(2,2);
sig=eye(2,2);
ar=armac(a,b,d,2,1,sig)

7.70   armax armax identification

[arc,la,lb,sig,resid]=armax(r,s,y,u,[b0f,prf])

• y : output process y(ny,n); ( ny: dimension of y , n : sample size)
• u : input process u(nu,n); ( nu: dimension of u , n : sample size)
• r and s : auto-regression orders r >=0 et s >=-1
• b0f : optional parameter. Its default value is 0 and it means that the coefficient b0 must be identified. if bof=1 the b0 is supposed to be zero and is not identified
• prf : optional parameter for display control. If prf =1, the default value, a display of the identified Arma is given.
• arc : a Scilab arma object (see armac)
• la : is the list(a,a+eta,a-eta) ( la = a in dimension 1) ; where eta is the estimated standard deviation. , a=[Id,a1,a2,...,ar] where each ai is a matrix of size (ny,ny)
• lb : is the list(b,b+etb,b-etb) (lb =b in dimension 1) ; where etb is the estimated standard deviation. b=[b0,.....,b_s] where each bi is a matrix of size (nu,nu)
• sig : is the estimated standard deviation of the noise and resid=[ sig*e(t0),....] (

 
   A(z^-1)y= B(z^-1)u + sig*e(t)

A(z) = 1+a1*z+...+a_r*z^r; ( r=0 => A(z)=1)
B(z) = b0+b1*z+...+b_s z^s ( s=-1 => B(z)=0)

 
[arc,a,b,sig,resid]=armax(); // will gives an example in dimension 1

7.71   armax1 armax identification

[a,b,d,sig,resid]=armax1(r,s,q,y,u,[b0f])

• y : output signal
• u : input signal
• r,s,q : auto regression orders with r >=0, s >=-1.
• b0f : optional parameter. Its default value is 0 and it means that the coefficient b0 must be identified. if bof=1 the b0 is supposed to be zero and is not identified
• a : is the vector [1,a1,...,a_r]
• b : is the vector [b0,......,b_s]
• d : is the vector [1,d1,....,d_q]
• sig : resid=[ sig*echap(1),....,];

   A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t)
   e(t) is a 1-dimensional white noise with variance 1.
   A(z)= 1+a1*z+...+a_r*z^r; ( r=0 => A(z)=1)
   B(z)= b0+b1*z+...+b_s z^s ( s=-1 => B(z)=0)
   D(z)= 1+d1*z+...+d_q*z^q  ( q=0 => D(z)=1)

7.72   arsimul armax simulation

[z]=arsimul(a,b,d,sig,u,[up,yp,ep])
[z]=arsimul(ar,u,[up,yp,ep])

• ar : an armax process. See armac.
• a : is the matrix[Id,a1,...,a_r] of dimension (n,(r+1)*n)
• b : is the matrix[b0,......,b_s] of dimension (n,(s+1)*m)
• d : is the matrix[Id,d_1,......,d_t] of dimension (n,(t+1)*n)
• u : is a matrix (m,N), which gives the entry u(:,j)=u_j
• sig : is a (n,n) matrix e_{k}is an n-dimensional Gaussian process with variance I
• up, yp : optional parameter which describe the past. up=[ u_0,u_{-1},...,u_{s-1}]; yp=[ y_0,y_{-1},...,y_{r-1}]; ep=[ e_0,e_{-1},...,e_{r-1}]; if they are omitted, the past value are supposed to be zero
• z : z=[ y(1),....,y(N)]

7.73   narsimul armax simulation ( using rtitr)

[z]=narsimul(a,b,d,sig,u,[up,yp,ep])
[z]=narsimul(ar,u,[up,yp,ep])

7.74   noisegen noise generation

[]=noisegen(pas,Tmax,sig)

noisegen(0.5,30,1.0);
x=-5:0.01:35;
y=feval(x,Noise);
plot(x,y);

7.75   odedi test of ode

[]=odedi()

7.76   prbs_a pseudo random binary sequences generation

[u]=prbs_a(n,nc,[ids])

u=prbs_a(50,10);
plot2d2("onn",(1:50)',u',1,"151",' ',[0,-1.5,50,1.5]);

7.77   reglin Linear regression

[a,b,sig]=reglin(x,y)

// simulation of data for a(3,5) and b(3,1)
x=rand(5,100);
aa=testmatrix('magi',5);aa=aa(1:3,:);
bb=[9;10;11]
y=aa*x +bb*ones(1,100)+ 0.1*rand(3,100);
// identification 
[a,b,sig]=reglin(x,y);
maxi(abs(aa-a))
maxi(abs(bb-b))
// an other example : fitting a polynom 
f=1:100; x=[f.*f; f];
y= [ 2,3]*x+ 10*ones(f) + 0.1*rand(f);
[a,b]=reglin(x,y)