filter_design.cpp
Go to the documentation of this file.
#include <itpp/signal/filter_design.h>
#include <itpp/signal/poly.h>
#include <itpp/signal/filter.h>
#include <itpp/signal/transforms.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/base/algebra/ls_solve.h>
#include <itpp/base/matfunc.h>
#include <itpp/base/specmat.h>
#include <itpp/base/math/trig_hyp.h>
#include <itpp/base/converters.h>
namespace itpp
{
void polystab(const vec &a, vec &out)
{
cvec r;
roots(a, r);
for (int i = 0; i < r.size(); i++) {
if (abs(r(i)) > 1)
r(i) = std::complex<double>(1.0) / conj(r(i));
}
out = real(std::complex<double>(a(0)) * poly(r));
}
void polystab(const cvec &a, cvec &out)
{
cvec r;
roots(a, r);
for (int i = 0; i < r.size(); i++) {
if (abs(r(i)) > 1)
r(i) = std::complex<double>(1.0) / conj(r(i));
}
out = a(0) * poly(r);
}
// ----------------------- freqz() ---------------------------------------------------------
void freqz(const cvec &b, const cvec &a, const int N, cvec &h, vec &w)
{
w = pi * linspace(0, N - 1, N) / double(N);
cvec ha, hb;
hb = fft(b, 2 * N);
hb = hb(0, N - 1);
ha = fft(a, 2 * N);
ha = ha(0, N - 1);
h = elem_div(hb, ha);
}
cvec freqz(const cvec &b, const cvec &a, const int N)
{
cvec h;
vec w;
freqz(b, a, N, h, w);
return h;
}
cvec freqz(const cvec &b, const cvec &a, const vec &w)
{
int la = a.size(), lb = b.size(), k = std::max(la, lb);
cvec h, ha, hb;
// Evaluate the nominator and denominator at the given frequencies
hb = polyval(zero_pad(b, k), to_cvec(cos(w), sin(w)));
ha = polyval(zero_pad(a, k), to_cvec(cos(w), sin(w)));
h = elem_div(hb, ha);
return h;
}
void freqz(const vec &b, const vec &a, const int N, cvec &h, vec &w)
{
w = pi * linspace(0, N - 1, N) / double(N);
cvec ha, hb;
hb = fft_real(b, 2 * N);
hb = hb(0, N - 1);
ha = fft_real(a, 2 * N);
ha = ha(0, N - 1);
h = elem_div(hb, ha);
}
cvec freqz(const vec &b, const vec &a, const int N)
{
cvec h;
vec w;
freqz(b, a, N, h, w);
return h;
}
cvec freqz(const vec &b, const vec &a, const vec &w)
{
int la = a.size(), lb = b.size(), k = std::max(la, lb);
cvec h, ha, hb;
// Evaluate the nominator and denominator at the given frequencies
hb = polyval(zero_pad(b, k), to_cvec(cos(w), sin(w)));
ha = polyval(zero_pad(a, k), to_cvec(cos(w), sin(w)));
h = elem_div(hb, ha);
return h;
}
void filter_design_autocorrelation(const int N, const vec &f, const vec &m, vec &R)
{
it_assert(f.size() == m.size(), "filter_design_autocorrelation: size of f and m vectors does not agree");
int N_f = f.size();
it_assert(f(0) == 0.0, "filter_design_autocorrelation: first frequency must be 0.0");
it_assert(f(N_f - 1) == 1.0, "filter_design_autocorrelation: last frequency must be 1.0");
// interpolate frequency-response
int N_fft = 512;
vec m_interp(N_fft + 1);
// unused variable:
// double df_interp = 1.0/double(N_fft);
m_interp(0) = m(0);
double inc;
int jstart = 0, jstop;
for (int i = 0; i < N_f - 1; i++) {
// calculate number of points to the next frequency
jstop = floor_i(f(i + 1) * (N_fft + 1)) - 1;
//std::cout << "jstart=" << jstart << "jstop=" << jstop << std::endl;
for (int j = jstart; j <= jstop; j++) {
inc = double(j - jstart) / double(jstop - jstart);
m_interp(j) = m(i) * (1 - inc) + m(i + 1) * inc;
}
jstart = jstop + 1;
}
vec S = sqr(concat(m_interp, reverse(m_interp(2, N_fft)))); // create a complete frequency response with also negative frequencies
R = ifft_real(to_cvec(S)); // calculate correlation
R = R.left(N);
}
// Calculate the AR coefficients of order \c n of the ARMA-process defined by the autocorrelation R
// using the deternined modified Yule-Walker method
// maxlag determines the size of the system to solve N>= n.
// If N>m then the system is overdetermined and a least squares solution is used.
// as a rule of thumb use N = 4*n
void modified_yule_walker(const int m, const int n, const int N, const vec &R, vec &a)
{
it_assert(m > 0, "modified_yule_walker: m must be > 0");
it_assert(n > 0, "modified_yule_walker: n must be > 0");
it_assert(N <= R.size(), "modified_yule_walker: autocorrelation function too short");
// create the modified Yule-Walker equations Rm * a = - rh
// see eq. (3.7.1) in Stoica and Moses, Introduction to spectral analysis
int M = N - m - 1;
mat Rm;
if (m - n + 1 < 0)
Rm = toeplitz(R(m, m + M - 1), reverse(concat(R(1, std::abs(m - n + 1)), R(0, m))));
else
Rm = toeplitz(R(m, m + M - 1), reverse(R(m - n + 1, m)));
vec rh = - R(m + 1, m + M);
// solve overdetermined system
a = backslash(Rm, rh);
// prepend a_0 = 1
a = concat(1.0, a);
// stabilize polynomial
a = polystab(a);
}
void arma_estimator(const int m, const int n, const vec &R, vec &b, vec &a)
{
it_assert(m > 0, "arma_estimator: m must be > 0");
it_assert(n > 0, "arma_estimator: n must be > 0");
it_assert(2*(m + n) <= R.size(), "arma_estimator: autocorrelation function too short");
// windowing the autocorrelation
int N = 2 * (m + n);
vec Rwindow = elem_mult(R.left(N), 0.54 + 0.46 * cos(pi * linspace(0.0, double(N - 1), N) / double(N - 1))); // Hamming windowing
// calculate the AR part using the overdetmined Yule-Walker equations
modified_yule_walker(m, n, N, Rwindow, a);
// --------------- Calculate MA part --------------------------------------
// use method in ref [2] section VII.
vec r_causal = Rwindow;
r_causal(0) *= 0.5;
vec h_inv_a = filter(1, a, concat(1.0, zeros(N - 1))); // see eq (50) of [2]
mat H_inv_a = toeplitz(h_inv_a, concat(1.0, zeros(m)));
vec b_causal = backslash(H_inv_a, r_causal);
// calculate the double-sided spectrum
int N_fft = 256;
vec H = 2.0 * real(elem_div(fft_real(b_causal, N_fft), fft_real(a, N_fft))); // calculate spectrum
// Do weighting and windowing in cepstrum domain
cvec cepstrum = log(to_cvec(H));
cvec q = ifft(cepstrum);
// keep only causal part of spectrum (windowing)
q.set_subvector(N_fft / 2, zeros_c(N_fft / 2));
q(0) *= 0.5;
cvec h = ifft(exp(fft(q))); // convert back to frequency domain, from cepstrum and do inverse transform to calculate impulse response
b = real(backslash(to_cmat(H_inv_a), h(0, N - 1))); // use Shank's method to calculate b coefficients
}
void yulewalk(const int N, const vec &f, const vec &m, vec &b, vec &a)
{
it_assert(f.size() == m.size(), "yulewalk: size of f and m vectors does not agree");
int N_f = f.size();
it_assert(f(0) == 0.0, "yulewalk: first frequency must be 0.0");
it_assert(f(N_f - 1) == 1.0, "yulewalk: last frequency must be 1.0");
vec R;
filter_design_autocorrelation(4*N, f, m, R);
arma_estimator(N, N, R, b, a);
}
} // namespace itpp
itpp::Array
General array class.
Definition array.h:105
itpp::Array::size
int size() const
Returns the number of data elements in the array object.
Definition array.h:155
converters.h
Definitions of converters between different vector and matrix types.
elem_math.h
Elementary mathematical functions - header file.
filter.h
Definitions of Filter classes and functions.
filter_design.h
Filter design functions.
it_assert
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
itpp::ifft
ITPP_EXPORT void ifft(const cvec &in, cvec &out)
Inverse Fast Fourier Transform.
itpp::fft_real
ITPP_EXPORT void fft_real(const vec &in, cvec &out)
Real Fast Fourier Transform.
itpp::ifft_real
ITPP_EXPORT void ifft_real(const cvec &in, vec &out)
Inverse Real Fast Fourier Transform.
itpp::fft
ITPP_EXPORT void fft(const cvec &in, cvec &out)
Fast Fourier Transform.
itpp::filter_design_autocorrelation
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)
Definition filter_design.cpp:153
itpp::yulewalk
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.
Definition filter_design.cpp:269
itpp::arma_estimator
void arma_estimator(const int m, const int n, const vec &R, vec &b, vec &a)
Estimation of ARMA model given the autocorrelation.
Definition filter_design.cpp:228
itpp::modified_yule_walker
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.
Definition filter_design.cpp:197
itpp::filter
vec filter(const vec &b, const vec &a, const vec &input)
ARMA filter function.
Definition filter.cpp:39
itpp::polystab
void polystab(const vec &a, vec &out)
Polynomial Stabilization.
Definition filter_design.cpp:45
itpp::freqz
void freqz(const cvec &b, const cvec &a, const int N, cvec &h, vec &w)
Frequency response of filter.
Definition filter_design.cpp:71
itpp::backslash
bool backslash(const mat &A, const vec &b, vec &x)
A general linear equation system solver.
Definition ls_solve.cpp:651
itpp::log
vec log(const vec &x)
The natural logarithm of the elements.
Definition log_exp.h:241
itpp::exp
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition log_exp.h:155
itpp::zero_pad
Vec< T > zero_pad(const Vec< T > &v, int n)
Zero-pad a vector to size n.
Definition matfunc.h:257
itpp::real
vec real(const cvec &data)
Real part of complex values.
Definition elem_math.cpp:157
itpp::sqr
vec sqr(const cvec &data)
Absolute square of elements.
Definition elem_math.cpp:36
itpp::conj
cvec conj(const cvec &x)
Conjugate of complex value.
Definition elem_math.h:226
itpp::polyval
vec polyval(const vec &p, const vec &x)
Evaluate polynomial.
Definition poly.cpp:135
itpp::poly
void poly(const vec &r, vec &p)
Create a polynomial of the given roots.
Definition poly.cpp:40
itpp::roots
void roots(const vec &p, cvec &r)
Calculate the roots of the polynomial.
Definition poly.cpp:66
itpp::reverse
Vec< T > reverse(const Vec< T > &in)
Reverse the input vector.
Definition matfunc.h:777
itpp::linspace
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition specmat.cpp:106
itpp::zeros
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
itpp::toeplitz
const cmat toeplitz(const cvec &c)
Generate symmetric Toeplitz matrix from vector c (complex valued)
Definition specmat.cpp:210
itpp::zeros_c
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
itpp::sin
vec sin(const vec &x)
Sine function.
Definition trig_hyp.h:54
itpp::cos
vec cos(const vec &x)
Cosine function.
Definition trig_hyp.h:58
ls_solve.h
Definitions of functions for solving linear equation systems.
matfunc.h
Various functions on vectors and matrices - header file.
itpp
itpp namespace
Definition itmex.h:37
itpp::to_cvec
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition converters.h:107
itpp::to_cmat
cmat to_cmat(const Mat< T > &m)
Converts a Mat<T> to cmat.
Definition converters.h:232
itpp::concat
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition array.h:486
itpp::pi
const double pi
Constant Pi.
Definition misc.h:103
itpp::elem_div
Mat< Num_T > elem_div(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise division of two matrices.
Definition mat.h:1688
itpp::floor_i
int floor_i(double x)
The nearest smaller integer.
Definition converters.h:350
itpp::elem_mult
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition mat.h:1582
itpp::abs
bin abs(const bin &inbin)
absolute value of bin
Definition binary.h:174
poly.h
Polynomial functions.
specmat.h
Definitions of special vectors and matrices.
transforms.h
Fourier, Hadamard, Walsh-Hadamard, and 2D Hadamard transforms - header file.
trig_hyp.h
Trigonometric and hyperbolic functions - header file.

Generated on Mon Apr 22 2024 09:45:27 for IT++ by Doxygen 1.9.8