Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 88 additions & 97 deletions lib/fir.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@

namespace dsplib {

//-------------------------------------------------------------------------------------------------
namespace {

template<class T>
static void _conv(const T* restrict x, const T* restrict h, T* restrict r, int nh, int nx) {
const int nr = nx - nh + 1;
Expand All @@ -21,6 +22,92 @@
}
}

arr_real _lowpass_fir(int n, real_t wn, const arr_real& win) {
if (win.size() != (n + 1)) {
DSPLIB_THROW("Window must be n+1 elements");

Check warning on line 27 in lib/fir.cpp

View check run for this annotation

Codecov / codecov/patch

lib/fir.cpp#L27

Added line #L27 was not covered by tests
}

const bool is_odd = (n % 2 == 1);
const int M = win.size();
const int L = M / 2;
const real_t fc = wn / 2;
const arr_real w = win.slice(0, L);

auto tt = arange(L) - real_t(n) / 2;
auto h = zeros(M);
h.slice(0, L) = sin(2 * pi * fc * tt) / tt * w;
if (!is_odd) {
h(L) = 2 * pi * fc;
h.slice(L + 1, M) = flip(h.slice(0, L));
} else {
h.slice(L, M) = flip(h.slice(0, L));
}

h /= sum(h);
return h;
}

arr_real _highpass_fir(int n, real_t wn, const arr_real& win) {
wn = 1 - wn;
int t1 = 0;
if (n % 2 == 1) {
n += 1;
t1 = 1;
}

auto h = _lowpass_fir(n, wn, win);
auto hh = arr_real(h.slice(t1, n, 2));
h.slice(t1, n, 2) = -hh;
return h;
}

arr_real _bandpass_fir(int n, real_t wn1, real_t wn2, const arr_real& win) {
wn1 = wn1 / 2;
wn2 = wn2 / 2;
auto wp = (wn2 - wn1) / 2;
auto wc = wn1 + wp;
auto h = _lowpass_fir(n, 2 * wp, win);
auto t = arange(h.size()) - real_t(n) / 2;
h = 2 * h * cos(2 * pi * wc * t);
return h;
}

arr_real _bandstop_fir(int n, real_t wn1, real_t wn2, const arr_real& win) {
if (n % 2 == 1) {
n += 1;
}

auto h = (-1) * _bandpass_fir(n, wn1, wn2, win);
h(n / 2) = h(n / 2) + 1;
return h;
}

bool _equal(const real_t& x1, const real_t& x2) {
return (abs(x1 - x2) < 2 * eps());
}

bool _is_symmetric(const dsplib::arr_real& h) {
const int n = h.size();
for (int i = 0; i < n / 2; i++) {
if (!_equal(h[i], h[n - i - 1])) {
return false;
}
}
return true;
}

bool _is_antisymmetric(const dsplib::arr_real& h) {
const int n = h.size();
for (int i = 0; i < n / 2; i++) {
if (!_equal(h[i], -h[n - i - 1])) {
return false;
}
}
return true;
}

} // namespace

//-------------------------------------------------------------------------------------------------
template<>
arr_real FirFilter<real_t>::conv(const arr_real& x, const arr_real& h) {
Expand All @@ -29,7 +116,6 @@
return r;
}

//-------------------------------------------------------------------------------------------------
template<>
arr_cmplx FirFilter<cmplx_t>::conv(const arr_cmplx& x, const arr_cmplx& h) {
arr_cmplx r(x.size() - h.size() + 1);
Expand All @@ -48,12 +134,10 @@
_x = complex(zeros(fft_len));
}

//-------------------------------------------------------------------------------------------------
FftFilter::FftFilter(const arr_real& h)
: FftFilter(complex(h)) {
}

//-------------------------------------------------------------------------------------------------
arr_cmplx FftFilter::process(const arr_cmplx& x) {
const int nr = (x.size() + _nx) / _n * _n;
arr_cmplx r(nr);
Expand All @@ -80,76 +164,11 @@
return r;
}

//-------------------------------------------------------------------------------------------------
arr_real FftFilter::process(const arr_real& x) {
//TODO: real optimization
return real(process(complex(x)));
}

//----------------------------------------------------------------------------------------------
static arr_real _lowpass_fir(int n, real_t wn, const arr_real& win) {
if (win.size() != (n + 1)) {
DSPLIB_THROW("Window must be n+1 elements");
}

const bool is_odd = (n % 2 == 1);
const int M = win.size();
const int L = M / 2;
const real_t fc = wn / 2;
const arr_real w = win.slice(0, L);

auto tt = arange(L) - real_t(n) / 2;
auto h = zeros(M);
h.slice(0, L) = sin(2 * pi * fc * tt) / tt * w;
if (!is_odd) {
h(L) = 2 * pi * fc;
h.slice(L + 1, M) = flip(h.slice(0, L));
} else {
h.slice(L, M) = flip(h.slice(0, L));
}

h /= sum(h);
return h;
}

//----------------------------------------------------------------------------------------------
static arr_real _highpass_fir(int n, real_t wn, const arr_real& win) {
wn = 1 - wn;
int t1 = 0;
if (n % 2 == 1) {
n += 1;
t1 = 1;
}

auto h = _lowpass_fir(n, wn, win);
auto hh = arr_real(h.slice(t1, n, 2));
h.slice(t1, n, 2) = -hh;
return h;
}

//----------------------------------------------------------------------------------------------
static arr_real _bandpass_fir(int n, real_t wn1, real_t wn2, const arr_real& win) {
wn1 = wn1 / 2;
wn2 = wn2 / 2;
auto wp = (wn2 - wn1) / 2;
auto wc = wn1 + wp;
auto h = _lowpass_fir(n, 2 * wp, win);
auto t = arange(h.size()) - real_t(n) / 2;
h = 2 * h * cos(2 * pi * wc * t);
return h;
}

//----------------------------------------------------------------------------------------------
static arr_real _bandstop_fir(int n, real_t wn1, real_t wn2, const arr_real& win) {
if (n % 2 == 1) {
n += 1;
}

auto h = (-1) * _bandpass_fir(n, wn1, wn2, win);
h(n / 2) = h(n / 2) + 1;
return h;
}

//----------------------------------------------------------------------------------------------
arr_real fir1(int n, real_t wn, FilterType ftype, const arr_real& win) {
assert(n > 0);
Expand All @@ -166,7 +185,6 @@
DSPLIB_THROW("Not supported for current filter type");
}

//----------------------------------------------------------------------------------------------
arr_real fir1(int n, real_t wn1, real_t wn2, FilterType ftype, const arr_real& win) {
assert(n > 0);
assert(wn1 < wn2);
Expand Down Expand Up @@ -200,33 +218,6 @@
return fir1(n, wn1, wn2, ftype, window::hamming(nn));
}

//----------------------------------------------------------------------------------------------
static bool _equal(const real_t& x1, const real_t& x2) {
return (abs(x1 - x2) < 2 * eps());
}

//----------------------------------------------------------------------------------------------
static bool _is_symmetric(const dsplib::arr_real& h) {
const int n = h.size();
for (int i = 0; i < n / 2; i++) {
if (!_equal(h[i], h[n - i - 1])) {
return false;
}
}
return true;
}

//----------------------------------------------------------------------------------------------
static bool _is_antisymmetric(const dsplib::arr_real& h) {
const int n = h.size();
for (int i = 0; i < n / 2; i++) {
if (!_equal(h[i], -h[n - i - 1])) {
return false;
}
}
return true;
}

//----------------------------------------------------------------------------------------------
FirType firtype(const dsplib::arr_real& h) {
const int n = h.size();
Expand Down
Loading