source: Sophya/trunk/ArchTOIPipe/Processors/dtoeplz.c@ 4047

Last change on this file since 4047 was 1944, checked in by aubourg, 24 years ago

decorrelateur wiener
File size: 1.3 KB
RevLine 
[1944]1#define NRANSI
2#include "nrutil.h"
3#define FREERETURN {free_dvector(h,1,n);free_dvector(g,1,n);return;}
4
5void dtoeplz(double r[], double x[], double y[], int n)
6{
7 int j,k,m,m1,m2;
8 double pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn;
9 double *g,*h;
10
11 if (r[n] == 0.0) nrerror("dtoeplz-1 singular principal minor");
12 g=dvector(1,n);
13 h=dvector(1,n);
14 x[1]=y[1]/r[n];
15 if (n == 1) FREERETURN
16 g[1]=r[n-1]/r[n];
17 h[1]=r[n+1]/r[n];
18 for (m=1;m<=n;m++) {
19 m1=m+1;
20 sxn = -y[m1];
21 sd = -r[n];
22 for (j=1;j<=m;j++) {
23 sxn += r[n+m1-j]*x[j];
24 sd += r[n+m1-j]*g[m-j+1];
25 }
26 if (sd == 0.0) nrerror("dtoeplz-2 singular principal minor");
27 x[m1]=sxn/sd;
28 for (j=1;j<=m;j++) x[j] -= x[m1]*g[m-j+1];
29 if (m1 == n) FREERETURN
30 sgn = -r[n-m1];
31 shn = -r[n+m1];
32 sgd = -r[n];
33 for (j=1;j<=m;j++) {
34 sgn += r[n+j-m1]*g[j];
35 shn += r[n+m1-j]*h[j];
36 sgd += r[n+j-m1]*h[m-j+1];
37 }
38 if (sd == 0.0 || sgd == 0.0) nrerror("dtoeplz-3 singular principal minor");
39 g[m1]=sgn/sgd;
40 h[m1]=shn/sd;
41 k=m;
42 m2=(m+1) >> 1;
43 pp=g[m1];
44 qq=h[m1];
45 for (j=1;j<=m2;j++) {
46 pt1=g[j];
47 pt2=g[k];
48 qt1=h[j];
49 qt2=h[k];
50 g[j]=pt1-pp*qt2;
51 g[k]=pt2-pp*qt1;
52 h[j]=qt1-qq*pt2;
53 h[k--]=qt2-qq*pt1;
54 }
55 }
56 nrerror("dtoeplz - should not arrive here!");
57}
58#undef FREERETURN
59#undef NRANSI
Note: See TracBrowser for help on using the repository browser.