← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to macros/signal/cspect.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/signal/cspect.sci
 
1
function [sm,cwp]=cspect(nlags,ntp,wtype,x,y,wpar)
 
2
//<sm,cwp>=cspect(nlags,ntp,wtype,x,y,wpar)
 
3
//Spectral estimation using the modified periodogram method.
 
4
//Cross-spectral estimate of x and y is calculated when both
 
5
//x and y are given.  Auto-spectral estimate of x is calculated
 
6
//if y is not given.
 
7
//  x        :data if vector, amount of input data if scalar
 
8
//  y        :data if vector, amount of input data if scalar
 
9
//  nlags    :number of correlation lags (pos. integer)
 
10
//  ntp      :number of transform points (pos. integer)
 
11
//  wtype    :window type ('re','tr','hm','hn','kr','ch')
 
12
//  wpar     :optional window parameters
 
13
//           :         for wtype='kr', wpar>0
 
14
//           :         for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
 
15
//  sm       :power spectral estimate in the interval [0,1]
 
16
//  cwp      :calculated value of unspecified Chebyshev
 
17
//           :window parameter
 
18
//
 
19
//!
 
20
//author: C. Bunks  date: 16 Sept 1988
 
21
// Copyright INRIA
 
22
 
 
23
   [lhs,rhs]=argn(0);
 
24
   cross=0;
 
25
 
 
26
//construct window
 
27
 
 
28
   if rhs==4 then,
 
29
      w=window(wtype,2*nlags-1);
 
30
   else if rhs==5 then,
 
31
      if wtype=='kr' then,
 
32
         wpar=y;
 
33
         w=window(wtype,2*nlags-1,wpar);
 
34
      else if wtype=='ch' then,
 
35
         wpar=y;
 
36
         [w,cwp]=window(wtype,2*nlags-1,wpar);
 
37
      else,
 
38
         cross=1;
 
39
         w=window(wtype,2*nlags-1);
 
40
      end,
 
41
      end,
 
42
   else,
 
43
      [w,cwp]=window(wtype,2*nlags-1,wpar);
 
44
      cross=1;
 
45
   end,
 
46
   end,
 
47
 
 
48
//estimate correlations
 
49
 
 
50
   if maxi(size(x))==1 then,
 
51
      nsects=int(x/(3*nlags));
 
52
      xlen=int(x/nsects);
 
53
      ss=0*ones(1,2*nlags);
 
54
      if cross==1 then,
 
55
         for k=1:nsects,
 
56
            xk=getx(xlen,1+(k-1)*xlen);
 
57
            yk=gety(xlen,1+(k-1)*xlen);
 
58
            ss=corr('update',xk,yk,ss);
 
59
         end,
 
60
         re=fft(ss,1)/x;
 
61
         re=[re(nlags:-1:1) re(2:nlags)];
 
62
      else,
 
63
         for k=1:nsects,
 
64
            xk=getx(xlen,1+(k-1)*xlen);
 
65
            ss=corr('update',xk,ss);
 
66
         end,
 
67
         re=fft(ss,1)/x;
 
68
         re=[re(nlags:-1:1) re(2:nlags)];
 
69
      end,
 
70
   else,
 
71
      if cross==1 then,
 
72
         [re,me]=corr(x,y,nlags);
 
73
         re=[re(nlags:-1:1) re(2:nlags)];
 
74
      else,
 
75
         [re,me]=corr(x,nlags);
 
76
         re=[re(nlags:-1:1) re(2:nlags)];
 
77
      end,
 
78
   end,
 
79
 
 
80
//window correlation estimate
 
81
 
 
82
   wre=w.*re;
 
83
 
 
84
//fourier transform to obtain spectral estimate
 
85
 
 
86
   wree=[wre 0*ones(1,ntp-2*nlags+1)];
 
87
   sm=abs(fft(wree,-1));
 
88
 
 
89
 
 
90