LocalCosine

Orthonormal local cosine bases (LocalCosineTransform), the cosine-packet tree with the Coifman-Wickerhauser best-basis search (CosinePackets) and the orthonormal DCT-IV they use (DCT4), declared in LocalCosine.hpp.

Conventions shared by the three classes:

  • Segment edges sit between samples: edge a separates sample a - 1 from sample a. A bell of half-width eps acts on the 2 eps samples a - eps .. a + eps - 1; 2 eps must not exceed the shortest segment, and eps = 0 is the block DCT-IV.

  • The rising cut-off is r(t) = sin(pi/4 (1 + beta(t))) at the sample centres t = (j + 1/2) / eps. The Coifman-Meyer bell takes beta as sin(pi t / 2) iterated three times; the sine bell takes beta(t) = t.

  • Each segment of M samples gets the orthonormal DCT-IV c_k = sqrt(2/M) sum_n y_n cos(pi/M (n + 1/2)(k + 1/2)), FFTW’s REDFT11 scaled by 1/sqrt(2M); it is its own inverse.

  • periodic window edges fold at 0 = N too, wrapping round, as the periodized wavelet transform does; free edges do not fold there. Both give orthonormal bases.

class LocalCosineTransform

A local cosine basis of segment length M on a window of N samples.

Forward writes the coefficients frequency-major: index k (N / M) + j holds DCT-IV bin k of segment j, so row k is the bin of centre (k + 1/2) fs / (2M), its N / M coefficients in time order, segment j covering samples [j M, (j + 1) M). That is the layout of a uniform wavelet-packet level of depth D = log2 M (2^D rows of N / 2^D), so the thresholds, the block rule and the trigger geometry of a packet level apply unchanged. ForwardSegments writes the same coefficients segment-major (index j M + k).

Public Types

enum Edge

What happens at the window edges 0 and N.

Values:

enumerator periodic

bells at the window edges too, wrapping round

enumerator free

no bell at the window edges

enum Bell

The shape of the rising cut-off r of every bell.

Values:

enumerator coifmanMeyer

r = sin(pi/4 (1 + s(s(s(t))))), s(t) = sin(pi t / 2)

enumerator sine

r = sin(pi/4 (1 + t))

Public Functions

LocalCosineTransform(unsigned int N, unsigned int M, int overlap = -1, enum Edge edge = periodic, enum Bell bell = coifmanMeyer)
Parameters:

  • N – window length, a multiple of M

  • M – segment length, a power of 2

  • overlap – bell half-width eps, at most M / 2; negative for M / 2

  • edge – what happens at the window edges

  • bell – the shape of the cut-off

Throws:

std::invalid_argument – on an inadmissible combination

void Forward(double *data)

In place, N values, frequency-major output.

void Inverse(double *data)

In place, N values, frequency-major input.

void ForwardSegments(double *data)

In place, N values, segment-major output.

void InverseSegments(double *data)

In place, N values, segment-major input.

inline unsigned int GetLength() const

The window length N.

inline unsigned int GetSegment() const

The segment length M.

inline unsigned int GetOverlap() const

The bell half-width eps.

inline enum Edge GetEdge() const

The treatment of the window edges.

inline enum Bell GetBell() const

The shape of the cut-off.

Public Static Functions

static std::vector<double> Ramp(unsigned int eps, enum Bell bell)

The cut-off r at the 2 eps sample centres t = (j + 1/2) / eps, j = -eps .. eps - 1.

static void Fold(double *data, unsigned int N, unsigned int a, const std::vector<double> &ramp, bool forward)

Fold (forward true) or unfold (forward false) the 2 eps samples about edge a of a window of N, in place; ramp is the output of Ramp(eps, bell). Indices wrap modulo N.

class CosinePackets

Cosine packets: the local cosine bases of every dyadic segmentation of the window, and the best basis among them (Coifman-Wickerhauser).

The tree has levels l = 0 .. levels - 1 of segment length maxSegment / 2^l, down to minSegment; node (l, i) is the segment [i m, (i + 1) m), m = maxSegment / 2^l, and its two children are its halves. Every node uses the same bell half-width eps at its edges, so an edge is folded the same way whatever level it belongs to: the coefficients of a node depend only on its own two edges, and any segmentation built from tree nodes (a “segmentation” below: segments in time order, each a tree node) is an orthonormal basis of the window.

The coefficients of a segmentation are written segment by segment in time order, each segment’s DCT-IV bins in frequency order.

The best basis minimises an additive cost, the sum over segments of a function of each segment’s coefficients, by the bottom-up search: a node is split when the best costs of its two children add up to strictly less than its own cost (a tie keeps the longer segment).

Public Types

enum Cost

The additive cost the best basis minimises.

Values:

enumerator l1

sum |c|

enumerator entropy

  • sum p ln p, p = c^2 / |x|^2 (Coifman-Wickerhauser)

enumerator wdfUniversal

  • sum (c/sigma)^2 over |c| > sqrt(2 ln N) sigma

enumerator wdfBlock

  • kept energy / sigma^2 of the block rule along each segment’s bins

Public Functions

CosinePackets(unsigned int N, unsigned int minSegment, unsigned int maxSegment = 0, int overlap = -1, enum LocalCosineTransform::Edge edge = LocalCosineTransform::periodic, enum LocalCosineTransform::Bell bell = LocalCosineTransform::coifmanMeyer)
Parameters:

  • N – window length, a power of 2

  • minSegment – shortest segment, a power of 2

  • maxSegment – longest segment, a power of 2 from minSegment to N; 0 for N

  • overlap – bell half-width, at most minSegment / 2; negative for minSegment / 2

  • edge – what happens at the window edges

  • bell – the shape of the cut-off

Throws:

std::invalid_argument – on an inadmissible combination

double BestBasis(const double *x, enum Cost cost, double sigma = 0.0)

The best basis of the window x (N values) for the cost. For the WDF costs sigma is the noise scale; zero or negative takes median |c| / 0.6745 over the coefficients of the finest level, which every basis of the tree shares for white noise. The block rule uses L = round(ln N) and lambda = 4.505, as WaveletThreshold; set them with SetBlock. l1 and entropy ignore sigma.

Returns:

the cost of the best basis

inline const std::vector<unsigned int> &GetSegmentation() const

The chosen segmentation: segment lengths in time order.

inline const std::vector<double> &GetCoefficients() const

The coefficients of x on the chosen segmentation (N values).

inline double GetSigma() const

The sigma the last BestBasis used (0 for l1 and entropy).

double GetNodeCost(unsigned int level, unsigned int index) const

Cost of node (level, index), from the last BestBasis.

Throws:

std::out_of_range – for a node outside the tree, or before the first BestBasis

double GetBestCost(unsigned int level, unsigned int index) const

Best cost of the subtree under node (level, index), from the last BestBasis.

Throws:

std::out_of_range – as GetNodeCost

double SegmentationCost(const std::vector<unsigned int> &segmentation) const

Sum of the node costs of a segmentation, from the last BestBasis.

Throws:

  • std::invalid_argument – when the segmentation is not made of tree nodes covering the window

  • std::out_of_range – before the first BestBasis

void Analyse(const double *x, const std::vector<unsigned int> &segmentation, double *out)

The coefficients of x (N values) on a segmentation, into out (N values): segment by segment in time order, each segment’s DCT-IV bins in frequency order.

Throws:

std::invalid_argument – when the segmentation is not made of tree nodes covering the window

void Inverse(const double *coeff, const std::vector<unsigned int> &segmentation, double *x)

The window (N values) from its coefficients on a segmentation, the exact inverse of Analyse.

Throws:

std::invalid_argument – as Analyse

double SegmentCost(const double *c, unsigned int M, enum Cost cost) const

The additive cost of the M coefficients of one segment, with the sigma and the window energy of the last BestBasis.

inline void SetBlock(unsigned int length, double lambda)

The parameters of the wdfBlock cost: the block length L in coefficients, 0 (the default) for round(ln N), and the energy threshold lambda in units of L sigma^2, 4.505 by default.

unsigned int GetBlockLength() const

The block length L in force.

inline double GetBlockLambda() const

The block threshold lambda in force.

inline unsigned int GetLength() const

The window length N.

inline unsigned int GetLevels() const

The number of levels of the tree, log2(maxSegment / minSegment) + 1.

unsigned int GetSegmentLength(unsigned int level) const

Segment length of a level, maxSegment / 2^level.

Throws:

std::out_of_range – for a level outside the tree

inline unsigned int GetOverlap() const

The bell half-width eps shared by every edge.

class DCT4

The orthonormal DCT-IV of a fixed length M, FFTW_REDFT11 scaled by 1/sqrt(2M). It is its own inverse.

Public Functions

explicit DCT4(unsigned int M)
Parameters:

M – transform length, positive

Throws:

  • std::invalid_argument – when M is 0

  • std::runtime_error – when FFTW cannot plan the transform

DCT4(const DCT4 &from)

A copy plans its own transform on its own buffer.

DCT4 &operator=(const DCT4 &from)

Assignment: the target takes the length and plans its own transform.

~DCT4()
void operator()(const double *in, double *out)

The transform of M values; in and out may be the same array.

inline unsigned int GetSize() const

The transform length M.