source: SH_SHM/trunk/source/polfiltr.c @ 344

Revision 16, 17.8 KB checked in by marcus, 15 years ago (diff)

r1 | svn | 2007-12-13 11:10:29 +0100 (Do, 13 Dez 2007) | 2 lines

Initial import
Line 
1
2/* file POLFILTR.C
3 *      ==========
4 *
5 * version 8, 22-May-2006
6 *
7 * polarisation filter of seismhandler
8 * K. Stammler, 7-APR-91
9 */
10
11
12/*
13 *
14 *  SeismicHandler, seismic analysis software
15 *  Copyright (C) 1992,  Klaus Stammler, Federal Institute for Geosciences
16 *                                       and Natural Resources (BGR), Germany
17 *
18 *  This program is free software; you can redistribute it and/or modify
19 *  it under the terms of the GNU General Public License as published by
20 *  the Free Software Foundation; either version 2 of the License, or
21 *  (at your option) any later version.
22 *
23 *  This program is distributed in the hope that it will be useful,
24 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
25 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
26 *  GNU General Public License for more details.
27 *
28 *  You should have received a copy of the GNU General Public License
29 *  along with this program; if not, write to the Free Software
30 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
31 *
32 */
33
34
35#include <stdio.h>
36#include <string.h>
37#include <math.h>
38#include "basecnst.h"
39#include BC_SYSBASE
40#include "shconst.h"
41#include "mxusrdef.h"
42
43
44#define COVMTXLTH 6
45        /* length of covariance matrix array */
46#define XX 0
47#define YY 1
48#define ZZ 2
49#define XY 3
50#define YZ 4
51#define XZ 5
52   /* index names of covariance matrix */
53
54typedef REAL COVMTX[COVMTXLTH];
55        /* covariance matrix */
56
57/* prototypes of local routines */
58static void pf_addcm( COVMTX a, COVMTX b );
59static void pf_subcm( COVMTX a, COVMTX b );
60static void pf_clrcm( COVMTX a );
61static void pf_buildcm( REAL x, REAL y, REAL z, COVMTX a );
62static void pf_evalcm( COVMTX c, REAL *l1, REAL *l2, REAL *l3 );
63static void pf_calccm( REAL x[], REAL y[], REAL z[], long li,
64        long hi, COVMTX c );
65static void pf_diagcm( COVMTX c, REAL l[MXC_DIM],
66        REAL trafo[MXC_DIM][MXC_DIM], int *big, STATUS *status );
67
68
69
70/*-----------------------------------------------------------------------*/
71
72
73
74void pf_filter( REAL xi[], REAL yi[], REAL zi[], REAL xo[],
75        REAL yo[], REAL zo[], long lth, int n, REAL p, int cmreset )
76
77/* filters 3-dim seismogram (x,y,z)
78 *
79 * parameters of routine
80 * REAL       xi[], yi[], zi[];   input; 3-dim input seismogram
81 * REAL       xo[], yo[], zo[];   output; filtered 3-dim seismogram
82 * long       lth;                input; length of traces
83 * int        n;                  input; coherence length
84 * REAL       p;                  input; exponent of filter factor
85 * int        cmreset;            input; the covariance matrix is
86 *                                       recomputed after cmreset steps
87 */
88{
89        /* local variables */
90        long     li, hi, mi; /* sample counter */
91        COVMTX   cm;         /* covariance matrix */
92        COVMTX   lm;         /* current matrix */
93        REAL     l1, l2, l3; /* eigenvalues */
94        REAL     fac;        /* filter factor */
95
96        /* executable code */
97
98        pf_clrcm( cm );
99        for  (li=(-n),hi=0; li<lth; li++,hi++)  {
100                if  (li >= 0)  {
101                        pf_buildcm( xi[li], yi[li], zi[li], lm );
102                        pf_subcm( cm, lm );
103                } /*endif*/
104                if  (hi < lth)  {
105                        pf_buildcm( xi[hi], yi[hi], zi[hi], lm );
106                        pf_addcm( cm, lm );
107                } /*endif*/
108                mi = (li+hi)/2;
109                if  (mi >=0  &&  mi < lth)  {
110                        if  ((mi+1) % cmreset == 0)      /* recompute cm to avoid */
111                                pf_calccm( xi, yi, zi, li>(-1) ? li+1 : 0, /* roundoff */
112                                        hi<lth ? hi : lth-1, cm );
113                        pf_evalcm( cm, &l1, &l2, &l3 );
114                        if  (l1 > SHC_EPSILON)  {
115                                fac = 1. - l2/l1;
116                                if  (p != 1.)  fac = pow(fac,p);
117                                xo[mi] = xi[mi]*fac;
118                                yo[mi] = yi[mi]*fac;
119                                zo[mi] = zi[mi]*fac;
120                        } else {
121                                xo[mi] = 0.;
122                                yo[mi] = 0.;
123                                zo[mi] = 0.;
124                        } /*endif*/
125                } /*endfor*/
126                sy_sharecpu();
127        } /*endfor*/
128
129} /* end of pf_filter */
130
131
132
133/*-----------------------------------------------------------------------*/
134
135
136
137void pf_vecfilter( REAL xi[], REAL yi[], REAL zi[], REAL xo[],
138        REAL yo[], REAL zo[], long lth, int n, REAL p, int cmreset, STATUS *status )
139
140/* filters 3-dim seismogram (x,y,z)
141 *
142 * parameters of routine
143 * REAL       xi[], yi[], zi[];   input; 3-dim input seismogram
144 * REAL       xo[], yo[], zo[];   output; filtered 3-dim seismogram
145 * long       lth;                input; length of traces
146 * int        n;                  input; coherence length
147 * REAL       p;                  input; exponent of filter factor
148 * int        cmreset;            input; the covariance matrix is
149 *                                       recomputed after cmreset steps
150 * STATUS     *status;            output; return status
151 */
152{
153        /* local variables */
154        long     li, hi, mi; /* sample counter */
155        COVMTX   cm;         /* covariance matrix */
156        COVMTX   lm;         /* current matrix */
157        REAL     fac;        /* filter factor */
158        REAL     l[MXC_DIM]; /* eigenvalues */
159        REAL     evec[MXC_DIM][MXC_DIM];   /* eigenvectors */
160        int      big;        /* index of biggest eigenvalue */
161        REAL     nvec[MXC_DIM]; /* normalized eigenvector of biggest value */
162        REAL     norm;       /* normalization constant */
163
164        /* executable code */
165
166        pf_clrcm( cm );
167        for  (li=(-n),hi=0; li<lth; li++,hi++)  {
168                if  (li >= 0)  {
169                        pf_buildcm( xi[li], yi[li], zi[li], lm );
170                        pf_subcm( cm, lm );
171                } /*endif*/
172                if  (hi < lth)  {
173                        pf_buildcm( xi[hi], yi[hi], zi[hi], lm );
174                        pf_addcm( cm, lm );
175                } /*endif*/
176                mi = (li+hi)/2;
177                if  (mi >=0  &&  mi < lth)  {
178                        if  ((mi+1) % cmreset == 0)      /* recompute cm to avoid */
179                                pf_calccm( xi, yi, zi, li>(-1) ? li+1 : 0, /* roundoff */
180                                        hi<lth ? hi : lth-1, cm );
181                        pf_diagcm( cm, l, evec, &big, status );
182                        if  (Severe(status))  return;
183                        if  (l[big] > SHC_EPSILON)  {
184                                /* eigenvector is evec[i][big], i=0,1,2 */
185                                /* normalize eigenvector */
186                                norm = sqrt( evec[0][big]*evec[0][big] +
187                                        evec[1][big]*evec[1][big] + evec[2][big]*evec[2][big] );
188                                if  (norm > SHC_EPSILON)  {
189                                        nvec[0] = fabs( evec[0][big] / norm );
190                                        nvec[1] = fabs( evec[1][big] / norm );
191                                        nvec[2] = fabs( evec[2][big] / norm );
192                                } else {
193                                        nvec[0] = nvec[1] = nvec[2] = 0.0;
194                                } /*endif*/
195                                if  (p != 1.)  {
196                                        nvec[0] = pow( nvec[0], p );
197                                        nvec[1] = pow( nvec[1], p );
198                                        nvec[2] = pow( nvec[2], p );
199                                } /*endif*/
200                                xo[mi] = xi[mi]*nvec[0];
201                                yo[mi] = yi[mi]*nvec[1];
202                                zo[mi] = zi[mi]*nvec[2];
203                        } else {
204                                xo[mi] = 0.;
205                                yo[mi] = 0.;
206                                zo[mi] = 0.;
207                        } /*endif*/
208                } /*endfor*/
209                sy_sharecpu();
210        } /*endfor*/
211
212} /* end of pf_vecfilter */
213
214
215
216/*-----------------------------------------------------------------------*/
217
218
219
220void pf_eigenvector( REAL xi[], REAL yi[], REAL zi[], REAL xo[],
221        REAL yo[], REAL zo[], long lth, int n, REAL p, int cmreset, STATUS *status )
222
223/* Returns normalized and weighted largest eigenvector for each sample position
224 *
225 * parameters of routine
226 * REAL       xi[], yi[], zi[];   input; 3-dim input seismogram
227 * REAL       xo[], yo[], zo[];   output; largest eigenvector
228 * long       lth;                input; length of traces
229 * int        n;                  input; coherence length
230 * REAL       p;                  input; exponent of filter factor
231 * int        cmreset;            input; the covariance matrix is
232 *                                       recomputed after cmreset steps
233 * STATUS     *status;            output; return status
234 */
235{
236        /* local variables */
237        long     li, hi, mi; /* sample counter */
238        COVMTX   cm;         /* covariance matrix */
239        COVMTX   lm;         /* current matrix */
240        REAL     fac;        /* filter factor */
241        REAL     l[MXC_DIM]; /* eigenvalues */
242        REAL     evec[MXC_DIM][MXC_DIM];   /* eigenvectors */
243        int      big;        /* index of biggest eigenvalue */
244        int      idx2, idx3; /* index of 2nd and 3rd largest eigenvector */
245        REAL     nvec[MXC_DIM]; /* normalized eigenvector of biggest value */
246        REAL     norm;       /* normalization constant */
247
248        /* executable code */
249
250        pf_clrcm( cm );
251        for  (li=(-n),hi=0; li<lth; li++,hi++)  {
252                if  (li >= 0)  {
253                        pf_buildcm( xi[li], yi[li], zi[li], lm );
254                        pf_subcm( cm, lm );
255                } /*endif*/
256                if  (hi < lth)  {
257                        pf_buildcm( xi[hi], yi[hi], zi[hi], lm );
258                        pf_addcm( cm, lm );
259                } /*endif*/
260                mi = (li+hi)/2;
261                if  (mi >=0  &&  mi < lth)  {
262                        if  ((mi+1) % cmreset == 0)      /* recompute cm to avoid */
263                                pf_calccm( xi, yi, zi, li>(-1) ? li+1 : 0, /* roundoff */
264                                        hi<lth ? hi : lth-1, cm );
265                        pf_diagcm( cm, l, evec, &big, status );
266                        if  (Severe(status))  return;
267                        switch  (big)  {
268                        case 0:  idx2 = 1;  idx3 = 2;  break;
269                        case 1:  idx2 = 0;  idx3 = 2;  break;
270                        default: idx2 = 0;  idx3 = 1;  break;
271                        } /*endswitch*/
272                        if  (fabs(l[idx2]) < fabs(l[idx3]))  idx2 = idx3;
273                        if  (fabs(l[big]) > SHC_EPSILON)  {
274                                /* eigenvector is evec[i][big], i=0,1,2 */
275                                /* normalize eigenvector */
276                                norm = sqrt( evec[0][big]*evec[0][big] +
277                                        evec[1][big]*evec[1][big] + evec[2][big]*evec[2][big] );
278                                if  (norm > SHC_EPSILON)  {
279                                        nvec[0] = fabs( evec[0][big] / norm );
280                                        nvec[1] = fabs( evec[1][big] / norm );
281                                        nvec[2] = fabs( evec[2][big] / norm );
282                                } else {
283                                        nvec[0] = nvec[1] = nvec[2] = 0.0;
284                                } /*endif*/
285                                /* weight eigenvecor by polarization */
286                                norm = 1.0 - l[idx2]/l[big];
287                                if  (p != 1.)
288                                        norm = pow( norm, p );
289                                xo[mi] = norm * nvec[0];
290                                yo[mi] = norm * nvec[1];
291                                zo[mi] = norm * nvec[2];
292                        } else {
293                                xo[mi] = 0.;
294                                yo[mi] = 0.;
295                                zo[mi] = 0.;
296                        } /*endif*/
297                } /*endfor*/
298                sy_sharecpu();
299        } /*endfor*/
300
301} /* end of pf_eigenvector */
302
303
304
305/*-----------------------------------------------------------------------*/
306
307
308
309void pf_poldir( REAL xi[], REAL yi[], REAL zi[], REAL xo[],
310        REAL yo[], REAL zo[], long lth, SAMPLE poldir[], REAL p )
311
312/* computes scalar product of each data point of the 3-dim input trace
313 * and a given vector poldir (both vectors are normalized before).  The
314 * data point then is weighted with the result powered by 'p'.
315 *
316 * parameters of routine
317 * REAL       xi[], yi[], zi[];   input; 3-dim input seismogram
318 * REAL       xo[], yo[], zo[];   output; filtered 3-dim seismogram
319 * long       lth;                input; length of traces
320 * int        poldir[];           input; polarization direction
321 * REAL       p;                  input; exponent of filter factor
322 */
323{
324        /* local variables */
325        REAL     norm;                 /* normalization constant */
326        SAMPLE   pol_z, pol_n, pol_e;  /* normalized polarization vector */
327        SAMPLE   cur_z, cur_n, cur_e;  /* normalized observed vector */
328        long     i;                    /* sample counter */
329
330        /* executable code */
331
332        /* normalize input vector */
333        norm = sqrt( poldir[0]*poldir[0] + poldir[1]*poldir[1]
334                + poldir[2]*poldir[2] );
335        if  (norm > SHC_EPSILON)  {
336                pol_z = poldir[0] / norm;
337                pol_n = poldir[1] / norm;
338                pol_e = poldir[2] / norm;
339        } else {
340                pol_z = pol_n = pol_e = 0.0;
341        } /*endif*/
342
343        for  (i=0; i<lth; i++)  {
344
345                norm = sqrt( xi[i]*xi[i] + yi[i]*yi[i] + zi[i]*zi[i] );
346                if  (norm > SHC_EPSILON)  {
347                        cur_z = xi[i] / norm;
348                        cur_n = yi[i] / norm;
349                        cur_e = zi[i] / norm;
350                } else {
351                        cur_z = cur_n = cur_e = 0.0;
352                } /*endif*/
353
354                /* compute powered scalar product */
355                norm = fabs( cur_z*pol_z + cur_n*pol_n + cur_e*pol_e );
356                if  (p != 1.0)  norm = pow( norm, p );
357
358                xo[i] = xi[i] * norm;
359                yo[i] = yi[i] * norm;
360                zo[i] = zi[i] * norm;
361
362        } /*endfor*/
363
364} /* end of pf_poldir */
365
366
367
368/*-----------------------------------------------------------------------*/
369
370
371
372void pf_bazdir( REAL xi[], REAL yi[], REAL zi[], REAL xo[],
373        REAL yo[], REAL zo[], long lth, REAL baz, REAL p )
374
375/* computes scalar product of each data point of the N,E input trace
376 * and a given back-azimuth (both vectors are normalized before).  The
377 * 3-dim data point then is weighted with the result powered by 'p'.
378 *
379 * parameters of routine
380 * REAL       xi[], yi[], zi[];   input; 3-dim input seismogram
381 * REAL       xo[], yo[], zo[];   output; filtered 3-dim seismogram
382 * long       lth;                input; length of traces
383 * REAL       baz;                input; back-azimuth in deg
384 * REAL       p;                  input; exponent of filter factor
385 */
386{
387        /* local variables */
388        REAL     norm;                 /* normalization constant */
389        SAMPLE   pol_n, pol_e;         /* normalized polarization vector */
390        SAMPLE   cur_n, cur_e;         /* normalized observed vector */
391        long     i;                    /* sample counter */
392
393        /* executable code */
394
395        /* get 2-dim input vector */
396        pol_n = cos( baz/SHC_RAD_TO_DEG );
397        pol_e = sin( baz/SHC_RAD_TO_DEG );
398
399        for  (i=0; i<lth; i++)  {
400
401                norm = sqrt( yi[i]*yi[i] + zi[i]*zi[i] );
402                if  (norm > SHC_EPSILON)  {
403                        cur_n = yi[i] / norm;
404                        cur_e = zi[i] / norm;
405                } else {
406                        cur_n = cur_e = 0.0;
407                } /*endif*/
408
409                /* compute powered scalar product */
410                norm = fabs( cur_n*pol_n + cur_e*pol_e );
411                if  (p != 1.0)  norm = pow( norm, p );
412
413                xo[i] = xi[i] * norm;
414                yo[i] = yi[i] * norm;
415                zo[i] = zi[i] * norm;
416
417        } /*endfor*/
418
419} /* end of pf_bazdir */
420
421
422
423/*-----------------------------------------------------------------------*/
424
425
426
427static void pf_addcm( COVMTX a, COVMTX b )
428
429/* adds matrix b to a
430 *
431 * parameters of routine
432 * COVMTX     a;       modify; b is added to this matrix
433 * COVMTX     b;       input; to be added
434 */
435{
436        /* local variables */
437        int      m;        /* element counter */
438
439        /* executable code */
440
441        for  (m=0; m<COVMTXLTH; m++)
442                a[m] += b[m];
443
444} /* end of pf_addcm */
445
446
447
448/*-----------------------------------------------------------------------*/
449
450
451
452static void pf_subcm( COVMTX a, COVMTX b )
453
454/* subtracts matrix b from a
455 *
456 * parameters of routine
457 * COVMTX     a;       modify; b is subtracted from this matrix
458 * COVMTX     b;       input; to be subtracted
459 */
460{
461        /* local variables */
462        int      m;        /* element counter */
463
464        /* executable code */
465
466        for  (m=0; m<COVMTXLTH; m++)
467                a[m] -= b[m];
468
469} /* end of pf_subcm */
470
471
472
473/*-----------------------------------------------------------------------*/
474
475
476
477static void pf_clrcm( COVMTX a )
478
479/* resets matrix a
480 *
481 * parameters of routine
482 * COVMTX     a;       output; to be reset
483 */
484{
485        /* local variables */
486        int      m;        /* element counter */
487
488        /* executable code */
489
490        for  (m=0; m<COVMTXLTH; m++)
491                a[m] = 0.;
492
493} /* end of pf_clrcm */
494
495
496
497/*-----------------------------------------------------------------------*/
498
499
500
501static void pf_buildcm( REAL x, REAL y, REAL z, COVMTX a )
502
503/* computes covariance matrix from 3-dim point coordinates
504 *
505 * parameters of routine
506 * REAL       x, y, z; input; coordinates of point
507 * COVMTX     a;       output; to be reset
508 */
509{
510        /* executable code */
511
512        a[XX] = x*x;
513        a[XY] = x*y;   a[YY] = y*y;
514        a[XZ] = x*z;   a[YZ] = y*z;   a[ZZ] = z*z;
515
516} /* end of pf_buildcm */
517
518
519
520/*-----------------------------------------------------------------------*/
521
522
523
524static void pf_calccm( REAL x[], REAL y[], REAL z[], long li,
525        long hi, COVMTX c )
526
527/* computes covariance matrix in sample window li..hi (inclusive
528 * both sample bounds)
529 *
530 * parameters of routine
531 * REAL       x[], y[], z[];    input; seismogram
532 * long       li, hi;           input; sample window
533 * COVMTX     c;                output; covariance matrix
534 */
535{
536        /* local variables */
537        long     i;           /* sample counter */
538        COVMTX   lc;          /* dyade */
539
540        /* executable code */
541
542        pf_clrcm( c );
543        for  (i=li; i<=hi; i++)  {
544                pf_buildcm( x[i], y[i], z[i], lc );
545                pf_addcm( c, lc );
546        } /*endif*/
547
548} /* end of pf_calccm */
549
550
551
552/*-----------------------------------------------------------------------*/
553
554
555static void pf_evalcm( COVMTX c, REAL *l1, REAL *l2, REAL *l3 )
556
557/* computes eigenvalues of matrix a
558 *
559 * parameters of routine
560 * COVMTX     c;             input; covariance matrix to be diagonalized
561 * REAL       *l1, *l2, *l3; output; eigenvalues
562 */
563{
564        /* local variables */
565        REAL     cp0, cp1, cp2;  /* coefficients of characteristic polynomial */
566        REAL     p, q;           /* corefficients of reduced equation */
567        REAL     D;              /* discriminante */
568        REAL     rho, phi;       /* scratch */
569
570        /* executable code */
571
572        /* characteristic polynomial: p(l) = l^3 + cp2*l^2 + cp1*l + cp0 */
573        cp2 = -c[XX] - c[YY] - c[ZZ];
574        cp1 = c[XX]*c[YY] + c[XX]*c[ZZ] + c[YY]*c[ZZ] -
575                        c[XY]*c[XY] - c[XZ]*c[XZ] - c[YZ]*c[YZ];
576        cp0 = c[XX]*c[YZ]*c[YZ] + c[YY]*c[XZ]*c[XZ] + c[ZZ]*c[XY]*c[XY] -
577                        c[XX]*c[YY]*c[ZZ] - 2.0*c[XY]*c[YZ]*c[XZ];
578
579        cp2 /= 3.;
580
581        /* reduced equation: l^3 + p*l  + q = 0 */
582        p = cp1 - 3.*cp2*cp2;
583        q = 2.*cp2*cp2*cp2 - cp2*cp1 + cp0;
584        if  (p > 0.)  {p = -p; printf( "*** this cannot happen ***\n" );}
585
586        rho = sqrt( -p );
587        rho = rho*rho*rho/5.1961524;    /* sqrt(27.) */
588        D = (rho > SHC_EPSILON) ? -q/(2.*rho) : 0.;
589        if  (D > 1.)  D = 1.;
590        if  (D < -1.)  D = -1.;
591        phi = acos( D );
592
593        /* compute eigenvalues */
594        D = (rho > 0.) ? 2.*pow(rho,1./3.) : 0.;
595        phi /= 3.;
596        *l1 = D * cos(phi);
597        *l2 = D * cos(phi+2.*SHC_PI/3.);
598        *l3 = D * cos(phi+4.*SHC_PI/3.);
599
600        /* undo reduction */
601        *l1 -= cp2;
602        *l2 -= cp2;
603        *l3 -= cp2;
604
605        /* sort ev's by size */
606        if  (*l1 < *l2)  {D = *l2; *l2 = *l1; *l1 = D;}
607        if  (*l2 < *l3)  {D = *l3; *l3 = *l2; *l2 = D;}
608        if  (*l1 < *l2)  {D = *l2; *l2 = *l1; *l1 = D;}
609
610} /* end of pf_evalcm */
611
612
613
614/*-----------------------------------------------------------------------*/
615
616
617
618static void pf_diagcm( COVMTX c, REAL l[MXC_DIM],
619        REAL trafo[MXC_DIM][MXC_DIM], int *big, STATUS *status )
620
621/* Returns eiganvalues and eigenvectors of symmetrical matrix c.
622 *
623 * parameters of routine
624 * COVMTX     c;                        input; covariance matrix
625 * REAL       l[MXC_DIM];               output; eigenvalues
626 * REAL       trafo[MXC_DIM][MXC_DIM];  output; eigenvectors
627 * int        *big;                     output; index of biggest value
628 * STATUS     *status;                  output; return status
629 */
630{
631        /* local variables */
632        REAL     sym_m[MXC_DIM][MXC_DIM];   /* symmetric input matrix */
633
634        /* executable code */
635
636        /* create input matrix */
637        sym_m[0][0] = c[XX];  sym_m[0][1] = c[XY];  sym_m[0][2] = c[XZ];
638        sym_m[1][0] = c[XY];  sym_m[1][1] = c[YY];  sym_m[1][2] = c[YZ];
639        sym_m[2][0] = c[XZ];  sym_m[2][1] = c[YZ];  sym_m[2][2] = c[ZZ];
640
641        /* diagonalize */
642        mx_diag_matrix( sym_m, l, trafo, status );
643
644        /* find biggest eigenvalue */
645        *big = 0;
646        if  (l[1] > l[*big])  *big = 1;
647        if  (l[2] > l[*big])  *big = 2;
648
649} /* end of pf_diagcm */
650
651
652
653/*-----------------------------------------------------------------------*/
Note: See TracBrowser for help on using the repository browser.