50 for (
int i = 0; i < r.size(); i++) {
52 r(i) = std::complex<double>(1.0) /
conj(r(i));
54 out =
real(std::complex<double>(a(0)) *
poly(r));
62 for (
int i = 0; i < r.size(); i++) {
64 r(i) = std::complex<double>(1.0) /
conj(r(i));
71 void freqz(
const cvec &b,
const cvec &a,
const int N, cvec &h, vec &w)
85 cvec
freqz(
const cvec &b,
const cvec &a,
const int N)
96 cvec
freqz(
const cvec &b,
const cvec &a,
const vec &w)
98 int la = a.size(), lb = b.size(), k =
std::max(la, lb);
111 void freqz(
const vec &b,
const vec &a,
const int N, cvec &h, vec &w)
125 cvec
freqz(
const vec &b,
const vec &a,
const int N)
130 freqz(b, a, N, h, w);
136 cvec
freqz(
const vec &b,
const vec &a,
const vec &w)
138 int la = a.size(), lb = b.size(), k =
std::max(la, lb);
155 it_assert(f.size() == m.size(),
"filter_design_autocorrelation: size of f and m vectors does not agree");
158 it_assert(f(0) == 0.0,
"filter_design_autocorrelation: first frequency must be 0.0");
159 it_assert(f(N_f - 1) == 1.0,
"filter_design_autocorrelation: last frequency must be 1.0");
163 vec m_interp(N_fft + 1);
170 int jstart = 0, jstop;
172 for (
int i = 0; i < N_f - 1; i++) {
174 jstop =
floor_i(f(i + 1) * (N_fft + 1)) - 1;
177 for (
int j = jstart; j <= jstop; j++) {
178 inc = double(j - jstart) / double(jstop - jstart);
179 m_interp(j) = m(i) * (1 - inc) + m(i + 1) * inc;
199 it_assert(m > 0,
"modified_yule_walker: m must be > 0");
200 it_assert(n > 0,
"modified_yule_walker: n must be > 0");
201 it_assert(N <= R.size(),
"modified_yule_walker: autocorrelation function too short");
214 vec rh = - R(m + 1, m + M);
230 it_assert(m > 0,
"arma_estimator: m must be > 0");
231 it_assert(n > 0,
"arma_estimator: n must be > 0");
232 it_assert(2*(m + n) <= R.size(),
"arma_estimator: autocorrelation function too short");
237 vec Rwindow =
elem_mult(R.left(N), 0.54 + 0.46 *
cos(
pi *
linspace(0.0,
double(N - 1), N) /
double(N - 1)));
244 vec r_causal = Rwindow;
250 vec b_causal =
backslash(H_inv_a, r_causal);
258 cvec q =
ifft(cepstrum);
261 q.set_subvector(N_fft / 2,
zeros_c(N_fft / 2));
269 void yulewalk(
const int N,
const vec &f,
const vec &m, vec &b, vec &a)
271 it_assert(f.size() == m.size(),
"yulewalk: size of f and m vectors does not agree");
274 it_assert(f(0) == 0.0,
"yulewalk: first frequency must be 0.0");
275 it_assert(f(N_f - 1) == 1.0,
"yulewalk: last frequency must be 1.0");
Mat< Num_T > elem_div(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise division of two matrices.
Various functions on vectors and matrices - header file.
void yulewalk(const int N, const vec &f, const vec &m, vec &b, vec &a)
ARMA filter design using a least-squares fit to the specified frequency-response. ...
Vec< T > reverse(const Vec< T > &in)
Reverse the input vector.
Definitions of Filter classes and functions.
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
cvec conj(const cvec &x)
Conjugate of complex value.
ITPP_EXPORT void ifft(const cvec &in, cvec &out)
Inverse Fast Fourier Transform.
cmat to_cmat(const Mat< T > &m)
Converts a Mat<T> to cmat.
ITPP_EXPORT void fft_real(const vec &in, cvec &out)
Real Fast Fourier Transform.
ITPP_EXPORT void fft(const cvec &in, cvec &out)
Fast Fourier Transform.
ITPP_EXPORT void ifft_real(const cvec &in, vec &out)
Inverse Real Fast Fourier Transform.
#define it_assert(t, s)
Abort if t is not true.
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
const cmat toeplitz(const cvec &c)
Generate symmetric Toeplitz matrix from vector c (complex valued)
vec sin(const vec &x)
Sine function.
void arma_estimator(const int m, const int n, const vec &R, vec &b, vec &a)
Estimation of ARMA model given the autocorrelation.
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
Definitions of converters between different vector and matrix types.
vec log(const vec &x)
The natural logarithm of the elements.
const double pi
Constant Pi.
vec exp(const vec &x)
Exp of the elements of a vector x.
T max(const Vec< T > &v)
Maximum value of vector.
Trigonometric and hyperbolic functions - header file.
void filter_design_autocorrelation(const int N, const vec &f, const vec &m, vec &R)
Calculate autocorrelation from the specified frequency-response (suitable for filter design) ...
vec filter(const vec &b, const vec &a, const vec &input)
ARMA filter functionThese functions implements a autoregressive moving average (ARMA) filter accordin...
void modified_yule_walker(const int m, const int n, const int N, const vec &R, vec &a)
Estimation of AR-part in an ARMA model given the autocorrelation.
vec polyval(const vec &p, const vec &x)
Evaluate polynomialEvaluate the polynomial p (of length at the points x The output is given by ...
bool backslash(const mat &A, const vec &b, vec &x)
A general linear equation system solver.
int floor_i(double x)
The nearest smaller integer.
void freqz(const cvec &b, const cvec &a, const int N, cvec &h, vec &w)
Frequency response of filter.
Definitions of functions for solving linear equation systems.
vec sqr(const cvec &data)
Absolute square of elements.
Definitions of special vectors and matrices.
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
void roots(const vec &p, cvec &r)
Calculate the roots of the polynomialCalculate the roots r of the polynomial p.
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
bin abs(const bin &inbin)
absolute value of bin
void polystab(const vec &a, vec &out)
Polynomial Stabilization.
vec cos(const vec &x)
Cosine function.
Elementary mathematical functions - header file.
vec real(const cvec &data)
Real part of complex values.
void poly(const vec &r, vec &p)
Create a polynomial of the given rootsCreate a polynomial p with roots r.
Vec< T > zero_pad(const Vec< T > &v, int n)
Zero-pad a vector to size n.
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.