/home/pclz-02/oschulz/Data/Science/Projects/Gerda/Software/blitzwave/src/Wavelet.h

// Copyright (C) 2003-2008 Oliver Schulz <oliver.schulz@tu-dortmund.de>

// 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.


#ifndef WAVELET_H
#define WAVELET_H

#include <cassert>
#include <cmath>
#include <limits>
#include <vector>
#include <string>
#include <iostream>
#include <blitz/array.h>

#include "arrayTools.h"

namespace bwave {



class Wavelet {
public:
        static bool equals(double x, double y);

public:

        class WaveletComponent {
        public:
                virtual ~WaveletComponent() { }
        };


        class LiftingStep : public WaveletComponent {
        public:
                enum Type { PRIMAL, DUAL, SSCALE, DSCALE };

        protected:
                Type m_type;
                int m_origin;
                std::vector<double> m_coeffs;
                double m_divisor;

                template<class tp_Type, int tp_size> inline void applyFS
                        (blitz::Array<tp_Type, 1> &s, blitz::Array<tp_Type, 1> &d, bool inverse, ExtensionMode be) const;

        public:
                Type type() const { return m_type; }

                int origin() const { return m_origin; }

                int size() const { return int(m_coeffs.size()); }

                double coeff(int i) const { return m_coeffs[i]; }

                double divisor() const { return m_divisor; }
                
                template<class tp_Type> inline void apply
                        (blitz::Array<tp_Type, 1> &s, blitz::Array<tp_Type, 1> &d, bool inverse, ExtensionMode be=CONSTANT_EXT) const;

                LiftingStep () {}

                LiftingStep (Type t, int o, double d, double c0);
                LiftingStep (Type t, int o, double d, double c0, double c1);
                LiftingStep (Type t, int o, double d, double c0, double c1, double c2);
                LiftingStep (Type t, int o, double d, double c0, double c1, double c2, double c3);
                LiftingStep (Type t, int o, double d, double c0, double c1, double c2, double c3, double c4);
                LiftingStep (Type t, int o, double d, double c0, double c1, double c2, double c3, double c4, double c5);
                
                bool operator==(const LiftingStep &s) const;
                
                virtual ~LiftingStep() {}
        };

protected:
        std::string m_spec;
        double m_nS, m_nD;
        std::vector<LiftingStep> m_steps;

public:
        const std::string& name() const { return m_spec; }

        bool integerCompatible() const;

        bool normalized() const;

        int steps() const { return int(m_steps.size()); }

        const LiftingStep& liftingStep(int i) const { return m_steps[i]; }

        double normS() const { return m_nS; }

        double normD() const { return m_nD; }

        bool operator==(const Wavelet &w) const;

        template<class tp_Type> void lift( blitz::Array<tp_Type, 1> &s,
                blitz::Array<tp_Type, 1> &d, ExtensionMode be=CONSTANT_EXT) const;

        template<class tp_Type> void invLift( blitz::Array<tp_Type, 1> &s,
                blitz::Array<tp_Type, 1> &d, ExtensionMode be=CONSTANT_EXT) const;

                
        blitz::Array<double, 1> forwardFkt(unsigned int size = 1<<9, int scale = 6,
                unsigned int trans = 3) const;
                                                                  
        blitz::Array<double, 1> inverseFkt(unsigned int size = 1<<9, int scale = 6,
                unsigned int trans = 3) const;
                                                                  
        Wavelet() {}

        Wavelet(const std::string &name, double ns, double nd, const LiftingStep &s0, const LiftingStep &s1);
        Wavelet(const std::string &name, double ns, double nd, const LiftingStep &s0, const LiftingStep &s1, const LiftingStep &s2);
        Wavelet(const std::string &name, double ns, double nd, const LiftingStep &s0, const LiftingStep &s1, const LiftingStep &s2, const LiftingStep &s3);

        virtual ~Wavelet() {}
};


std::ostream & operator<<(std::ostream &os, const Wavelet::LiftingStep &s);

std::ostream & operator<<(std::ostream &os, const Wavelet &w);


extern const Wavelet WL_CDF_1_1;

extern const Wavelet WL_CDF_2_2;

extern const Wavelet WL_CDF_3_1;

extern const Wavelet WL_CDF_3_3;

extern const Wavelet WL_CDF_4_2;

extern const Wavelet WL_CDF_97;

extern const Wavelet WL_HAAR;

extern const Wavelet WL_LEG_5_3;

extern const Wavelet WL_CUBIC_SPLINE;

extern const Wavelet WL_D_4;


template<class tp_Type, int tp_size> void Wavelet::LiftingStep::applyFS
        (blitz::Array<tp_Type, 1> &s, blitz::Array<tp_Type, 1> &d, bool inverse, ExtensionMode be) const
{
        assert(size() == tp_size);
        tp_Type coeffs[tp_size];
        for (int i=0; i<size(); ++i) coeffs[i] = (tp_Type) coeff(i);
        GenFilter<tp_Type, tp_size> liftFilter(origin(), coeffs, (tp_Type) divisor());
                
        if (type() == PRIMAL) {
                if (!inverse) liftFilter.applyAdd(d, s, be);
                else liftFilter.applySub(d, s, be);
        } else if (type() == DUAL) {
                if (!inverse) liftFilter.applyAdd(s, d, be);
                else liftFilter.applySub(s, d, be);
        } else if (type() == SSCALE) {
                if (!inverse) {
                        s *= tp_Type(coeff(0));
                        if (divisor()!=1.) s /= tp_Type(divisor()); 
                }
                else {
                        if (divisor()!=1.) s *= tp_Type(divisor()); 
                        s /= tp_Type(coeff(0));
                }
        } else if (type() == DSCALE) {
                if (!inverse) {
                        d *= tp_Type(coeff(0));
                        if (divisor()!=1.) d /= tp_Type(divisor()); 
                }
                else {
                        if (divisor()!=1.) d *= tp_Type(divisor()); 
                        d /= tp_Type(coeff(0));
                }
        } else {
                assert(false);
        }
}


template<class tp_Type> void Wavelet::LiftingStep::apply
        (blitz::Array<tp_Type, 1> &s, blitz::Array<tp_Type, 1> &d, bool inverse, ExtensionMode be) const
{
        switch(size()) {
                case 1: applyFS<tp_Type, 1>(s, d, inverse, be); break;
                case 2: applyFS<tp_Type, 2>(s, d, inverse, be); break;
                case 3: applyFS<tp_Type, 3>(s, d, inverse, be); break;
                //case 4: applyFS<tp_Type, 4>(s, d, inverse, be); break;
                //case 5: applyFS<tp_Type, 5>(s, d, inverse, be); break;
                //case 6: applyFS<tp_Type, 6>(s, d, inverse, be); break;
                default: assert(false);
        }
}


template<class tp_Type> void Wavelet::lift( blitz::Array<tp_Type, 1> &s,
        blitz::Array<tp_Type, 1> &d, ExtensionMode be) const
{
        for (int step=0; step<steps(); ++step)
                liftingStep(step).apply(s, d, false, be);
}


template<class tp_Type> void Wavelet::invLift( blitz::Array<tp_Type, 1> &s,
        blitz::Array<tp_Type, 1> &d, ExtensionMode be) const
{
        for (int step=steps()-1; step>=0; --step)
                liftingStep(step).apply(s, d, true, be);
}


} // namespace bwave

#endif