cseife.c 4.68 KB
Newer Older
Simone Butzer's avatar
Simone Butzer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
/*! \file cseife.c
 * \brief make seife functions available within C code (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id: cseife.c 3381 2010-11-15 12:33:00Z tforb $
 * \author Thomas Forbriger
 * \date 14/01/2005
 * 
 * make seife functions available within C code (implementation)
 * 
 * Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 14/01/2005   V1.0   Thomas Forbriger
 *  - 11/07/2005   V1.1   support debug mode
 *  - 15/11/2010   V1.2   do not use tfmacros.h
 * 
 * ============================================================================
 */
#define TF_CSEIFE_C_VERSION \
  "TF_CSEIFE_C   V1.2"
#define TF_CSEIFE_C_CVSID \
  "$Id: cseife.c 3381 2010-11-15 12:33:00Z tforb $"

#include <cseife.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* initialize global flags */
struct seife_flags seife_global_flags={ 0 };

/* enter debug mode */
void seife_debug_mode_on()
{
  seife_global_flags.debug=1;
}

/* Butterworth lowpass (period t0, order o) */
void seife_lpb(double* x, int n, double dt, double t0, int o)
{
  seife_butterworth(x, n, dt, IIRlpb, t0, o);
}

/* Butterworth highpass (period t0, order o) */
void seife_hpb(double* x, int n, double dt, double t0, int o)
{
  seife_butterworth(x, n, dt, IIRhpb, t0, o);
}

/* 2nd order lowpass (period t0, damping h) */
void seife_lp2(double* x, int n, double dt, double t0, double h)
{
  seife_rekfl(x, x, n, seife_rfk(IIRlp2, t0/dt, h, 0., 0.));
}

/* 2nd order highpass (period t0, damping h) */
void seife_hp2(double* x, int n, double dt, double t0, double h)
{
  seife_first(x, n);
  seife_rekfl(x, x, n, seife_rfk(IIRhp2, t0/dt, h, 0., 0.));
}

/* 2nd order bandpass (period t0, damping h) */
void seife_bp2(double* x, int n, double dt, double t0, double h)
{
  seife_rekfl(x, x, n, seife_rfk(IIRbp2, t0/dt, h, 0., 0.));
}

/* 1st order lowpass (period t0) */
void seife_lp1(double* x, int n, double dt, double t0)
{
  seife_rekfl(x, x, n, seife_rfk(IIRlp1, t0/dt, 0., 0., 0.));
}

/* 1st order highpass (period t0) */
void seife_hp1(double* x, int n, double dt, double t0)
{
  seife_first(x, n);
  seife_rekfl(x, x, n, seife_rfk(IIRhp1, t0/dt, 0., 0., 0.));
}

/* integration (time constant t0) */
void seife_int(double* x, int n, double dt, double t0)
{
  double t00=t0;
  if (t00 < 1.e-22) { t00=1.; }
  seife_rekfl(x, x, n, seife_rfk(IIRint, t00/dt, 0., 0., 0.));
}

/* 1st order highpass equalizer (former period t0s, new period t0) */
void seife_he1(double* x, int n, double dt, double t0s, double t0)
{
  seife_rekfl(x, x, n, seife_rfk(IIRhe1, t0/dt, 0., t0s/dt, 0.));
}

/* 1st order lowpass equalizer (former period t0s, new period t0) */
void seife_le1(double* x, int n, double dt, double t0s, double t0)
{
  struct seife_IIRcoef coef=seife_rfk(IIRle1, t0/dt, 0., t0s/dt, 0.);
  double fac=t0s/t0;
  coef.f0 *= fac;
  coef.f1 *= fac;
  seife_rekfl(x, x, n, coef);
}

/* 2nd order highpass equalizer (former: period t0s and damping hs,
 *                               new: period t0 and damping h) 
 */
void seife_he2(double* x, int n, double dt, 
               double t0s, double hs, double t0, double h)
{
  if (seife_global_flags.debug)
  {
    fprintf(stderr, 
            "DEBUG (seife_he2): n=%d, dt=%f, t0s=%f, hs=%f, t0=%f, h=%f\n",
            n, dt, t0s, hs, t0, h);
  }
  seife_rekfl(x, x, n, seife_rfk(IIRhe2, t0/dt, h, t0s/dt, hs));
}

/* 2nd order lowpass equalizer (former: period t0s and damping hs,
 *                              new: period t0 and damping h) 
 */
void seife_le2(double* x, int n, double dt, 
               double t0s, double hs, double t0, double h)
{
  struct seife_IIRcoef coef=seife_rfk(IIRle2, t0/dt, h, t0s/dt, hs);
  double fac=t0s/t0;
  fac *= fac;
  coef.f0 *= fac;
  coef.f1 *= fac;
  coef.f2 *= fac;
  seife_rekfl(x, x, n, coef);
}

/* set baseline to first value */
void seife_first(double* x, int n)
{
  double x0=x[0];
  int i;
  for (i=0; i<n; ++i) { x[i] -= x0; }
}

/* general Butterworth filter framework */
void seife_butterworth(double* x, int n, double dt,
                       enum seife_filter_type IIRtype, 
                       double t0, int m)
{
  enum seife_filter_type t1st, t2nd;
  if (IIRtype == IIRhpb)
  {
    seife_first(x, n);
    t1st=IIRhp1;
    t2nd=IIRhp2;
  }
  else if (IIRtype == IIRlpb)
  {
    t1st=IIRlp1;
    t2nd=IIRlp2;
  }
  else
  {
    SEIFE_CHECKERROR( 1, "seife_butterworth", "illegal IIRtype" )
  }
  int mm=m/2;
  if (m > (2*mm))
  { seife_rekfl(x, x, n, seife_rfk(t1st, t0/dt, 0., 0., 0.)); }
  double pih=2.*atan(1.);
  int j;
  for (j=1; j<=mm; ++j)
  {
    double h=sin(pih*(2*j-1)/m);
    seife_rekfl(x, x, n, seife_rfk(t2nd, t0/dt, h, 0., 0.)); 
  }
}


/* ----- END OF cseife.c ----- */