function [sp,df] = tsmooth(p,kw) % % s = tsmooth(p,kw) computes a spectral estimate by smoothing periodogram % p, using a TRIANGULAR window of full (odd) width = kw points % % see pdgm.m for details of suitable periodogram calculation. % see also rsmooth.m, csmooth.m for box & cosine-bell smoothing windows % % MIM Feb/90 % May-92 a call [sp,df] = tsmooth(p,kw) now gives an extra o/p vector, df, % the first element of which is the no. deg. of freedom in the Chi-square % approx. to the dist. of the spectral estimates, sp. % if kw < 1, kw=1;msg='warning, zero length window call',end; kworig=kw; % % generate triangular window of full (odd) width = original % kw (if it wasn't odd go to next higher odd no.). window ->0 % at lags +/- (kw/2+1) kw=fix(kw/2);k0=kw+1; w=ones(2*kw+1,1); for k=1:kw, w(k0+k)=1-k/(kw+1); w(k0-k)=1-k/(kw+1); end; w=w/sum(w); %%msg='check window',keyboard; % % % % reflect periodogram about end points ..... % original periodogram now starts with zero freq. pt at index = kw+1 % Note use of conjugate in reflected bits so we can handle cross-periodograms % This won't affect auto-periodograms. np=length(p); pp=ones(np+(2*kw),1); pp(kw+1:kw+np)=p(1:np); pp(kw:-1:1)=conj(pp(kw+1+1:kw+1+kw)); pp(kw+np+1:kw+np+kw)=conj(pp(kw+np-1:-1:kw+np-kw)); pp(kw+1)=0;%should already be so if mean correction done at pdgm stage; %%msg='check reflection in pp',keyboard; % % % run window along pp... sp=conv(pp,w); % % select desired bit..... sp = sp(2*kw:2*kw+np-1); % % adjust first kw+1 points for zero in periodogram at f=0 for k=0:kw,fix=1/(1-w(k0+k));sp(k+1)=sp(k+1)*fix;end; df=(1+0*sp)*(2/sum(w.*w));