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
|
/*!
* \file
* \brief Filter design functions
* \author Tony Ottosson
*
* -------------------------------------------------------------------------
*
* IT++ - C++ library of mathematical, signal processing, speech processing,
* and communications classes and functions
*
* Copyright (C) 1995-2008 (see AUTHORS file for a list of contributors)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* -------------------------------------------------------------------------
*/
#ifndef FILTER_DESIGN_H
#define FILTER_DESIGN_H
#include <itpp/base/vec.h>
namespace itpp {
/*!
\addtogroup filters
*/
/*!
\brief Polynomial Stabilization
\ingroup filters
\author Tony Ottosson
Stabilizes the polynomial transfer function by replacing all roots outside
the unit cirlce with their reflection inside the unit circle.
@{
*/
void polystab(const vec &a, vec &out);
inline vec polystab(const vec &a) { vec temp; polystab(a, temp); return temp; }
void polystab(const cvec &a, cvec &out);
inline cvec polystab(const cvec &a) { cvec temp; polystab(a, temp); return temp; }
/*! @} */
/*!
\brief Frequency response of filter
\ingroup filters
\author Tony Ottosson
Calculates the N-point frequency response of the supplied digital filter over the frequencies w.
If w is not given the response is evaluated over the range 0 to \f$\pi\f$ with N values.
The default value of N is 512.
If \c w is supplied polyval() is used. Otherwise the calculation is based on the fft.
@{
*/
void freqz(const cvec &b, const cvec& a, const int N, cvec &h, vec &w);
cvec freqz(const cvec &b, const cvec& a, const int N = 512);
cvec freqz(const cvec &b, const cvec& a, const vec &w);
void freqz(const vec &b, const vec& a, const int N, cvec &h, vec &w);
cvec freqz(const vec &b, const vec& a, const int N = 512);
cvec freqz(const vec &b, const vec& a, const vec &w);
/*! @} */
/*!
\brief Calculate autocorrelation from the specified frequency-response (suitable for filter design)
\ingroup filters
\author Tony Ottosson
Calculates the autocorrelation function of size \c N corresponding to the specified frequency response.
Useful as a first step in designing filters.
The vectors \c f and \c m is the frequency response. The frequencies should be between 0 and 1.0
(equal to half the sample rate) in increasing order. Both 0.0 and 1.0 must be included. The frequency
response is upsampled to 512 points and the autocorrelation is ifft of the power
magnitude response of the upsampled frequency response.
*/
void filter_design_autocorrelation(const int N, const vec &f, const vec &m, vec &R);
/*!
\brief Estimation of AR-part in an ARMA model given the autocorrelation
\ingroup filters
\author Tony Ottosson
Estimates the AR-part of an ARMA model from the given autocorrelation. The AR part is of order \c n.
The overdetermined modified Yule-Walker equations are used.
If \f$N>n\f$ then the system is overdetermined and a least squares solution is used.
As a rule of thumb use \f$N = 4 n\f$
The parameter \c m is the order of the MA-part such that \f$R(k) = 0, \forall \|k\| > m\f$.
The supplied autocorrelation should at least be of size \c N.
References:
Stoica and Moses, Introduction to spectral analysis, Prentice Hall, 1997.
*/
void modified_yule_walker(const int m, const int n, const int N, const vec &R, vec &a);
/*!
\brief Estimation of ARMA model given the autocorrelation
\ingroup filters
\author Tony Ottosson
Estimates an ARMA model from the given autocorrelation. The AR part is of order \c n and the MA part
is of order \c m.
The AR part (the denominator) is calcuated using the modified Yule-Walker equations. The
the MA part (the nominator) is calculated by calculating the inverse magnitude spectrum using FFTs of
size 512 which is an AR-system. This AR-system is then solved using the Levinson-Durbin algorithm.
The supplied autocorrelation is windowed using a Hamming window of size \f$2(m+n)\f$ and hence should at
least be of that size.
References:
[1] Stoica and Moses, Introduction to spectral analysis, Prentice Hall, 1997.
[2] B. Friedlander and B. Porat, The modified Yule-Walker method of ARMA spectral estimation,
IEEE Trans. Aerospace and Electronic Systems, Vol. AES-20, No. 2, pp. 158--173, March 1984.
*/
void arma_estimator(const int m, const int n, const vec &R, vec &b, vec &a);
/*!
\brief ARMA filter design using a least-squares fit to the specified frequency-response
\ingroup filters
\author Tony Ottosson
The arma_estimator() function is used to calculate the a and b coefficients.
The vectors \c f and \c m is the frequency response. The frequencies should be between 0 and 1.0
(equal to half the sample rate) in increasing order. Both 0.0 and 1.0 must be included. The
filter_design_autocorrelation() fucnction is used to interpolate the frequency response and calculate
the corresponding autocorrelation.
Observe: this function will not always give exactly the same result as the matlab yulewalk function.
*/
void yulewalk(const int N, const vec &f, const vec &m, vec &b, vec &a);
} // namespace itpp
#endif // #ifndef FILTER_DESIGN_H
|
