shownew.m

to find stastical moments

function[spect2,f2]=shownew(b,a,n)

%function[spect2,f2]=shownew(b,a,n)

%  This function computes the "spectrum" (i.e., the squared

%  amplitude of the frequency response function) of a stochastic process

%

%       a(L)x(t) = b(L) e(t)

%

%  where e(t) is a white noise with unit variance and

%  a(L) = 1 +a1*L + ... + am*L^m, b(L) = b0 + b1*L + ... + br*L^r.

%  The frequency response of x to e is computed using FREQ.  Then

%  the spectrum is computed as the squared amplitude of the frequency

%  response.

%

%  The parameter n must be a positive integer power of 2.  The parameter

%  n determines the number of frequencies at which the spectrum is computed.

%  The parameter "st" is a string with which the program will label

%  the spectrum during plotting.

%

%  The spectrum and frequencies at which the spectrum is computed are

%  returned in spect and f, respectively.

s=freq(b,a,n);

spect=(abs(s)).^2;

f=2*pi*(0:n-1)/n;

n2=n/2;

f2=f(1:n2);

spect2=log(spect(1:n2));

%axis ([0 pi -1 1]);

plot(f2,spect2)

axis ([0 pi -inf inf])