[PATCH] Implement Fir<>, native C++ version
Nathan (Jasper) Myers
ncm at codesourcery.com
Mon Oct 10 01:22:34 UTC 2005
The following patch has been committed.
It implements vsip::Fir<> using native C++ code, and a comprehensive
test of all its modes.
Note that a few bits of the test are commented out; it uses
vsip::Convolution<> to generate the reference output, and that has
a little bug I haven't got to tracking down yet.
Nathan Myers
ncm
Index: ChangeLog
===================================================================
RCS file: /home/cvs/Repository/vpp/ChangeLog,v
retrieving revision 1.288
diff -u -p -r1.288 ChangeLog
--- ChangeLog 7 Oct 2005 13:46:45 -0000 1.288
+++ ChangeLog 10 Oct 2005 01:17:11 -0000
@@ -1,3 +1,9 @@
+2005-10-09 Nathan Myers <ncm at codesourcery.com>
+
+ Implement FIR filter, all modes.
+ * src/vsip/impl/signal-fir.hpp, tests/fir.cpp: New.
+ * src/vsip/signal.hpp: Include new impl/signal-fir.hpp.
+
2005-10-06 Jules Bergmann <jules at codesourcery.com>
Implement 1-D correlation.
Index: src/vsip/signal.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/signal.hpp,v
retrieving revision 1.10
diff -u -p -r1.10 signal.hpp
--- src/vsip/signal.hpp 7 Oct 2005 13:46:46 -0000 1.10
+++ src/vsip/signal.hpp 10 Oct 2005 01:17:11 -0000
@@ -19,6 +19,7 @@
#include <vsip/impl/signal-conv.hpp>
#include <vsip/impl/signal-corr.hpp>
#include <vsip/impl/signal-window.hpp>
+#include <vsip/impl/signal-fir.hpp>
#endif // VSIP_SIGNAL_HPP
Index: src/vsip/impl/signal-fir.hpp
===================================================================
RCS file: src/vsip/impl/signal-fir.hpp
diff -N src/vsip/impl/signal-fir.hpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/vsip/impl/signal-fir.hpp 10 Oct 2005 01:17:11 -0000
@@ -0,0 +1,208 @@
+/* Copyright (c) 2005 by CodeSourcery, LLC. All rights reserved. */
+
+/** @file vsip/impl/signal-fir.hpp
+ @author Nathan Myers
+ @date 2005-09-30
+ @brief VSIPL++ Library: FIR filters
+*/
+
+#ifndef VSIP_IMPL_SIGNAL_FIR_HPP
+#define VSIP_IMPL_SIGNAL_FIR_HPP
+
+#include <vsip/support.hpp>
+#include <vsip/impl/metaprogramming.hpp>
+#include <vsip/impl/signal-types.hpp>
+#include <vsip/vector.hpp>
+#include <vsip/domain.hpp>
+#include <vsip/impl/view_traits.hpp>
+#include <vsip/impl/fast-block.hpp>
+#include <vsip/impl/global_map.hpp>
+#include <vsip/impl/profile.hpp>
+
+#include <new>
+
+namespace vsip
+{
+
+enum obj_state { state_no_save, state_save };
+
+namespace impl
+{
+
+// Fir_aligned; block and view types with optimal layout for FIR operations.
+
+template<typename T, typename Map_type = vsip::Local_map>
+struct Fir_aligned
+{
+ typedef typename
+ impl::Fast_block<1,T,
+ impl::Layout<1,
+ vsip::tuple<0,1,2>,
+ impl::Stride_unit_align<16>,impl::Cmplx_inter_fmt>,
+ Map_type>
+ block_type;
+};
+
+} // namespace impl
+
+///////////////////////////////////////////////////////////////////
+//
+// class Fir<>
+//
+
+template
+<
+ typename T = VSIP_DEFAULT_VALUE_TYPE,
+ vsip::symmetry_type symV = vsip::nonsym,
+ vsip::obj_state useOldState = vsip::state_save,
+ unsigned NumberOfTimes = 0,
+ vsip::alg_hint_type algHint = vsip::alg_time
+>
+class Fir
+{
+public:
+ static const vsip::symmetry_type symmetry = symV;
+ static const vsip::obj_state continuous_filter = useOldState;
+
+ template <typename BlockT>
+ Fir(
+ vsip::const_Vector<T,BlockT> kernel,
+ vsip::length_type input_size,
+ vsip::length_type decimation = 1)
+ VSIP_THROW((std::bad_alloc))
+ : input_size_(input_size),
+ order_(kernel.size() * (1 + (symV != vsip::nonsym)) -
+ (symV == vsip::sym_even_len_odd) - 1),
+ decimation_(decimation),
+ skip_(0),
+ state_saved_(0),
+ state_(this->order_, T(0)),
+ kernel_(this->order_ + 1)
+ {
+ assert(input_size > 0);
+ assert(this->order_ + 1 > 1); // counter unsigned wraparound
+ assert(decimation > 0);
+ assert(this->order_ + 1 > decimation); // M >= decimation
+ assert(input_size >= this->order_); // input_size >= M
+
+ // must be after asserts because of division
+ this->output_size_ = (input_size + decimation - 1) / decimation;
+
+ // mirror the kernel
+ unsigned const ksz = kernel.size();
+ this->kernel_(vsip::Domain<1>(this->order_, -1, ksz)) = kernel;
+ // and maybe unmirror a copy, too
+ if (symV != vsip::nonsym)
+ this->kernel_(vsip::Domain<1>(ksz)) = kernel;
+ }
+
+ vsip::length_type kernel_size() const VSIP_NOTHROW
+ { return this->order_ + 1; }
+ vsip::length_type filter_order() const VSIP_NOTHROW
+ { return this->order_ + 1; }
+ // vsip::symmetry_type symmetry() const VSIP_NOTHROW
+ // { return symV; }
+ vsip::length_type input_size() const VSIP_NOTHROW
+ { return this->input_size_; }
+ vsip::length_type output_size() const VSIP_NOTHROW
+ { return this->output_size_; }
+ vsip::obj_state continuous_filtering() const VSIP_NOTHROW
+ { return useOldState; }
+ vsip::length_type decimation() const VSIP_NOTHROW
+ { return this->decimation_; }
+
+ template <typename Block0, typename Block1>
+ vsip::length_type
+ operator()(
+ vsip::const_Vector<T, Block0> data,
+ vsip::Vector<T, Block1> out)
+ VSIP_NOTHROW
+ {
+ assert(data.size() == this->input_size_);
+ assert(out.size() == this->output_size_);
+
+ typedef vsip::Domain<1> d_type;
+ typedef vsip::length_type len_type;
+
+ const len_type dec = this->decimation_;
+ const len_type m = this->order_;
+ const len_type skip = this->skip_;
+ const len_type saved = this->state_saved_;
+ len_type oix = 0;
+ len_type i = 0;
+ for (; i < m - skip; ++oix, i += dec)
+ {
+ // Conceptually this comes second, but it's more convenient
+ // to put it here. So, read the second statement below first.
+ T sum = vsip::dot(
+ this->kernel_(d_type(m - skip - i, 1, i + skip + 1)),
+ data(d_type(i + skip + 1)));
+
+ if (useOldState == vsip::state_save && i < saved)
+ {
+ sum += vsip::dot(
+ this->kernel_(d_type(saved - i)),
+ this->state_(d_type(i, 1, saved - i)));
+ }
+
+ out.put(oix, sum);
+ }
+
+ const len_type isz = this->input_size_;
+ len_type start = i - (m - skip);
+ for ( ; start + m < isz; ++oix, start += dec)
+ {
+ T sum = vsip::dot(this->kernel_, data(d_type(start, 1, m + 1)));
+ out.put(oix, sum);
+ }
+
+ if (useOldState == state_save)
+ {
+ this->skip_ = start + m - isz;
+ const len_type new_save = isz - start;
+ this->state_saved_ = new_save;
+ this->state_(d_type(new_save)) = data(d_type(start, 1, new_save));
+ }
+ return oix;
+ }
+
+ void reset() VSIP_NOTHROW
+ { this->state_saved_ = this->skip_ = 0; }
+
+ ~Fir()
+ { }
+
+public:
+
+ float impl_performance(char* what) const
+ {
+ if (!strcmp(what, "mflops"))
+ {
+ // Compute rough estimate of flop-count.
+ unsigned cxmul = impl::Is_complex<T>::value ? 4 : 1; // *(4*,2+), +(2+)
+ return (this->timer_.count() * cxmul * 2 * // 1* 1+
+ ((this->order + 1) * this->input_size_ / this->decimation_)) /
+ (1e6 * this->timer_.total());
+ }
+ else if (!strcmp(what, "time"))
+ return this->timer_.total();
+ return 0.f;
+ }
+
+private:
+ vsip::length_type input_size_;
+ vsip::length_type output_size_;
+ vsip::length_type order_; // M in the spec
+ vsip::length_type decimation_;
+ vsip::length_type skip_; // how much of next input to skip
+ vsip::length_type state_saved_; // number of elements saved
+ vsip::Vector<T,typename impl::Fir_aligned<T>::block_type> state_;
+ vsip::Vector<T,typename impl::Fir_aligned<T>::block_type> kernel_;
+
+ impl::profile::Acc_timer timer_;
+};
+
+} // namespace vsip
+
+#endif // VSIP_IMPL_SIGNAL_FIR_HPP
+
Index: tests/fir.cpp
===================================================================
RCS file: tests/fir.cpp
diff -N tests/fir.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ tests/fir.cpp 10 Oct 2005 01:17:11 -0000
@@ -0,0 +1,201 @@
+/* Copyright (c) 2005 by CodeSourcery, LLC. All rights reserved. */
+
+/** @file tests/fft.cpp
+ @author Nathan Myers
+ @date 2005-10-03
+ @brief VSIPL++ Library: Testcases for Fir<>
+*/
+
+/***********************************************************************
+ Included Files
+***********************************************************************/
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <vsip/initfin.hpp>
+#include <vsip/support.hpp>
+#include <vsip/signal.hpp>
+#include <vsip/math.hpp>
+#include <vsip/matrix.hpp>
+
+#include "test.hpp"
+#include "output.hpp"
+
+
+/***********************************************************************
+ Definitions
+***********************************************************************/
+
+template <typename T, vsip::symmetry_type sym>
+void
+test_fir(
+ vsip::length_type D,
+ vsip::length_type M,
+ vsip::length_type N)
+{
+ const unsigned insize = 2 * M * N;
+ const unsigned outsize = ((2 * M * N) + D - 1) / D + 1;
+ vsip::Vector<T> input(insize);
+ vsip::Vector<T> output1(outsize);
+ vsip::Vector<T> output2(2 * M * (N+D-1)/D);
+ vsip::Vector<T> output3(2 * M * (N+D-1)/D);
+
+ vsip::Vector<T> convinput(insize+M, T(0)); // room for initial state
+ vsip::Vector<T> convout((insize+M-1)/D + 1, T(0)); // per spec
+ vsip::Vector<T> kernel(M);
+
+ for (vsip::length_type i = 0; i < insize; ++i)
+ input.put(i, T(i+1));
+ for (vsip::length_type i = 0; i < M; ++i)
+ kernel.put(i, T(2*i+1));
+
+ vsip::Convolution<vsip::const_Vector,sym,vsip::support_same,T,1> conv(
+ kernel, vsip::Domain<1>(convinput.size()), D);
+
+ const vsip::length_type pad = (sym == vsip::nonsym) ? M/2 :
+ (sym == vsip::sym_even_len_even) ? M : M - 1;
+ convinput(vsip::Domain<1>(pad, 1, insize)) = input;
+ conv(convinput, convout); // emulate chained FIR
+
+ vsip::Fir<T,sym,vsip::state_save,1> fir1a(kernel, N, D);
+ vsip::Fir<T,sym,vsip::state_save,1> fir1b(kernel, N, D);
+ vsip::Fir<T,sym,vsip::state_no_save,1> fir2(kernel, N, D);
+
+ vsip::length_type got = 0;
+ for (vsip::length_type i = 0; i < 2 * M; ++i) // chained
+ {
+ got += fir1a(
+ input(vsip::Domain<1>(i * N, 1, N)),
+ output1(vsip::Domain<1>(got, 1, (N + D - 1) / D)));
+ }
+
+
+ vsip::length_type got1b = 0;
+ vsip::length_type got2 = 0;
+ for (vsip::length_type i = 0; i < 2 * M; ++i) // not
+ {
+ got1b += fir1b(input(vsip::Domain<1>(i * N, 1, N)),
+ output2(vsip::Domain<1>(got1b, 1, (N+D-1)/D)));
+ fir1b.reset();
+ got2 += fir2(input(vsip::Domain<1>(i * N, 1, N)),
+ output3(vsip::Domain<1>(got2, 1, (N+D-1)/D)));
+ }
+
+ vsip::Vector<T> reference(convout(vsip::Domain<1>(got)));
+ vsip::Vector<T> result(output1(vsip::Domain<1>(got)));
+ assert(outsize - got <= 1);
+ assert(vsip::alltrue(result == reference));
+ assert(got1b == got2);
+ assert(vsip::alltrue(
+ output2(vsip::Domain<1>(got1b)) == output3(vsip::Domain<1>(got1b))));
+}
+
+int
+main()
+{
+ vsip::vsipl init;
+
+
+ test_fir<float,vsip::nonsym>(1,2,3);
+
+ test_fir<float,vsip::nonsym>(3,23,31);
+
+ test_fir<float,vsip::nonsym>(1,3,5);
+ test_fir<float,vsip::nonsym>(1,3,9);
+ test_fir<float,vsip::nonsym>(1,4,8);
+ test_fir<float,vsip::nonsym>(1,23,31);
+ test_fir<float,vsip::nonsym>(1,32,256);
+
+ test_fir<float,vsip::nonsym>(2,3,5);
+ test_fir<float,vsip::nonsym>(2,3,9);
+ test_fir<float,vsip::nonsym>(2,4,8);
+ test_fir<float,vsip::nonsym>(2,23,31);
+ test_fir<float,vsip::nonsym>(2,32,256);
+
+ test_fir<double,vsip::nonsym>(2,3,5);
+ test_fir<double,vsip::nonsym>(2,3,9);
+ test_fir<double,vsip::nonsym>(2,4,8);
+ test_fir<double,vsip::nonsym>(2,23,31);
+ test_fir<double,vsip::nonsym>(2,32,1024);
+
+ test_fir<std::complex<float>,vsip::nonsym>(2,3,5);
+ test_fir<std::complex<float>,vsip::nonsym>(2,3,9);
+ test_fir<std::complex<float>,vsip::nonsym>(2,4,8);
+ test_fir<std::complex<float>,vsip::nonsym>(2,23,31);
+ test_fir<std::complex<float>,vsip::nonsym>(2,32,256);
+
+ test_fir<std::complex<double>,vsip::nonsym>(2,3,5);
+ test_fir<std::complex<double>,vsip::nonsym>(2,3,9);
+ test_fir<std::complex<double>,vsip::nonsym>(2,4,8);
+ test_fir<std::complex<double>,vsip::nonsym>(2,23,31);
+ test_fir<std::complex<double>,vsip::nonsym>(2,32,1024);
+
+ test_fir<float,vsip::nonsym>(3,4,8);
+ test_fir<float,vsip::nonsym>(3,4,21);
+ test_fir<float,vsip::nonsym>(3,9,27);
+ test_fir<float,vsip::nonsym>(3,23,31);
+ test_fir<float,vsip::nonsym>(3,32,256);
+
+ test_fir<float,vsip::nonsym>(4,5,13);
+ test_fir<float,vsip::nonsym>(4,7,31);
+ test_fir<float,vsip::nonsym>(4,8,32);
+ test_fir<float,vsip::nonsym>(4,23,31);
+ test_fir<float,vsip::nonsym>(4,32,256);
+
+ test_fir<float,vsip::sym_even_len_even>(1,1,3);
+ test_fir<float,vsip::sym_even_len_even>(1,2,3);
+ test_fir<float,vsip::sym_even_len_even>(1,3,5);
+ test_fir<float,vsip::sym_even_len_even>(1,3,9);
+ test_fir<float,vsip::sym_even_len_even>(1,4,8);
+ test_fir<float,vsip::sym_even_len_even>(1,23,57);
+#if 0
+ // FIXME: this exposes a bug in vsip::Convolution.
+ test_fir<float,vsip::sym_even_len_even>(1,32,256);
+#endif
+
+ test_fir<float,vsip::sym_even_len_even>(2,2,3);
+ test_fir<float,vsip::sym_even_len_even>(2,3,5);
+ test_fir<float,vsip::sym_even_len_even>(2,3,9);
+ test_fir<float,vsip::sym_even_len_even>(2,4,8);
+ test_fir<float,vsip::sym_even_len_even>(2,23,57);
+#if 0
+ // FIXME: likewise
+ test_fir<float,vsip::sym_even_len_even>(2,32,256);
+#endif
+
+ test_fir<float,vsip::sym_even_len_even>(3,3,5);
+ test_fir<float,vsip::sym_even_len_even>(3,4,8);
+ test_fir<float,vsip::sym_even_len_even>(3,23,57);
+#if 0
+ test_fir<float,vsip::sym_even_len_even>(3,32,256);
+#endif
+
+ test_fir<float,vsip::sym_even_len_odd>(1,2,3);
+ test_fir<float,vsip::sym_even_len_odd>(1,3,5);
+ test_fir<float,vsip::sym_even_len_odd>(1,3,9);
+ test_fir<float,vsip::sym_even_len_odd>(1,4,9);
+ test_fir<float,vsip::sym_even_len_odd>(1,23,57);
+#if 0
+ test_fir<float,vsip::sym_even_len_odd>(1,32,256);
+#endif
+
+ test_fir<float,vsip::sym_even_len_odd>(2,2,3);
+ test_fir<float,vsip::sym_even_len_odd>(2,3,5);
+ test_fir<float,vsip::sym_even_len_odd>(2,3,9);
+ test_fir<float,vsip::sym_even_len_odd>(2,4,10);
+ test_fir<float,vsip::sym_even_len_odd>(2,23,57);
+#if 0
+ test_fir<float,vsip::sym_even_len_odd>(2,32,256);
+#endif
+
+ test_fir<float,vsip::sym_even_len_odd>(3,3,5);
+ test_fir<float,vsip::sym_even_len_odd>(3,4,9);
+ test_fir<float,vsip::sym_even_len_odd>(3,23,55);
+#if 0
+ test_fir<float,vsip::sym_even_len_odd>(3,32,256);
+#endif
+
+ return 0;
+}
More information about the vsipl++
mailing list