mlpt

Multiscale local 1-D polynomial transform

collapse all in page

Syntax

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x,t)
[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x,t,numLevel)
[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x)
[coefs,T,coefsPerLevel,scalingMoments] = mlpt(___,Name=Value)

Description

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x,t) returns the multiscale local polynomial 1-D transform (MLPT) of input signal x sampled at the sampling instants, t. If x or t contain NaNs, the union of the NaNs in x and t is removed before obtaining the mlpt.

example

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x,t,numLevel) returns the transform for numLevel resolution levels.

example

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(x) uses uniform sampling instants for x as the time instants if x does not contain NaNs. If x contains NaNs, the NaNs are removed from x and the nonuniform sampling instants are obtained from the numeric elements of x.

example

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(___,Name=Value) specifies mlpt properties using one or more name-value arguments and any of the previous syntaxes. For example, mlpt(x,DualMoments=4) specifies four dual vanishing moments.

example

Examples

collapse all

Open Live Script

Create a signal with nonuniform sampling and verify good reconstruction.

Create and plot a sine wave with nonuniform sampling.

timeVector = 0:0.01:1;
sineWave = sin(2*pi*timeVector);

samplesToErase = randi(100,1,100);
sineWave(samplesToErase) = [];
timeVector(samplesToErase) = [];

plot(timeVector,sineWave,"o")
grid on
title("Signal")
xlabel("Time (s)")
ylabel("Amplitude")

Figure contains an axes object. The axes object with title Signal, xlabel Time (s), ylabel Amplitude contains a line object which displays its values using only markers.

Use mlpt to obtain the multiscale local 1-D polynomial transform of the signal. Visualize the coefficients.

[coefs,T,coefsPerLevel,scalingMoments] = mlpt(sineWave,timeVector);
stem(coefs)
grid on
title("Wavelet Coefficients")

Figure contains an axes object. The axes object with title Wavelet Coefficients contains an object of type stem.

Use imlpt to obtain the inverse multiscale local 1-D polynomial transform of the coefficients. Plot the original signal and the reconstruction.

y = imlpt(coefs,T,coefsPerLevel,scalingMoments);

plot(timeVector,sineWave,"o")
hold on
plot(T,y,"x")
hold off
grid on
legend("Original Signal","Reconstruction")
xlabel("Time (s)")
ylabel("Amplitude")

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Original Signal, Reconstruction.

Inspect the total error to verify good reconstruction.

reconstructionError = sum(abs(y'-sineWave))
reconstructionError = 
3.0951e-15
Open Live Script

Specify nondefault dual moments by using the mlpt function. Compare the results of analysis and synthesis using the default and nondefault dual moments.

Create an input signal and visualize it.

t = linspace(0,1,200)';
x = cos(10*pi*t.^2);
plot(t,x)
title("Signal")

Figure contains an axes object. The axes object with title Signal contains an object of type line.

By default, DualMoments, the number of dual vanishing moments, is 2. Perform the forward and inverse transform for the input signal twice, with DualMoments set to 2 and 3.

[w2,t2,nj2,scalingmoments2] = mlpt(x,t);
y2 = imlpt(w2,t2,nj2,scalingmoments2);

[w3,t3,nj3,scalingmoments3] = mlpt(x,t,DualMoments=3);
y3 = imlpt(w3,t3,nj3,scalingmoments3,DualMoments=3);

Plot both reconstructions.

plot(t,x)
hold on
plot(t2,y2,'o')
plot(t3,y3,'*')
hold off
legend("Original", ...
       "DualMoments = 2", ...
       "DualMoments = 3");

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent Original, DualMoments = 2, DualMoments = 3.

For each reconstruction, compute the mean error. Verify perfect reconstruction.

fprintf("Mean Error\nDualMoments = 2: %e\nDualMoments = 3: %e", ...
    mean(abs(y2-x)),mean(abs(y3-x)))
Mean Error
DualMoments = 2: 1.539914e-16
DualMoments = 3: 7.797582e-17
Open Live Script

Resolution levels are the number of cascaded local polynomial smoothing operations. The details at each resolution level are obtained by predicting one half the samples based on a local polynomial interpolation of the other half. The difference between the predicted and actual values are the details at each resolution level. The scaling coefficients at each coarser resolution level are smoother versions of the higher resolution scaling coefficients. Only the final-level scaling coefficients are retained.

Increasing the number of resolution levels enables you to analyze narrowband coefficients for a computational and memory cost.

Create a dual-tone input signal, x, that contains high and low frequencies.

fs = 1000;
t = (0:1/fs:10)';
x = sin(499*pi.*t) + sin(2*pi.*t);

Use mlpt to obtain coefficients for minimum and maximum resolution levels. Print the computation time.

tic
[w1,~,nj1,m1] = mlpt(x,t,1);
computationTime1 = toc;
disp("Level one computation time: "+computationTime1+" seconds")
Level one computation time: 0.93461 seconds

tic
[w13,~,nj13,m13] = mlpt(x,t,13);
computationTime13 = toc;
disp("Level thirteen computation time: "+computationTime13+" seconds")
Level thirteen computation time: 1.3304 seconds
Open Live Script

If your time instants are not known or specified, you can calculate the MLPT using default time instants.

Load a data signal corrupted with NaNs and with unknown time instants. Calculate the MLPT without specifying time instants. The resulting implied time instants is a vector of valid indices of the corrupted signal.

load CorruptedData

[w,t,nj,scalingMoments] = mlpt(yCorrupt);

Calculate the inverse MLPT and visualize the results. Reinsert NaNs to visualize gaps in the signal.

z = imlpt(w,t,nj,scalingMoments);

zToPlot = NaN(numel(yCorrupt),1);
zToPlot(t) = z;

plot(yCorrupt,LineWidth=2.5)
hold on
plot(zToPlot,LineWidth=1)
hold off
legend("Original Signal","Reconstructed Signal")
xlabel("Time Instants")

Figure contains an axes object. The axes object with xlabel Time Instants contains 2 objects of type line. These objects represent Original Signal, Reconstructed Signal.

Input Arguments

collapse all

Input signal, specified as a vector or matrix.

  • matrix — x must have at least two rows. mlpt operates independently on each column of x. The number of elements in t must equal the row dimension of x. Any NaNs in the columns of x must occur in the same rows.

  • vector — x and t must have the same number of elements.

Data Types: double

Sampling instants corresponding to the input signal, specified as a vector, duration array, or datetime array of monotonically increasing real values. The default value depends on the length of the input signal, x.

Data Types: double | duration | datetime

Number of resolution levels, specified as a positive integer. The maximum value of numLevel depends on the shape of the input signal, x:

  • matrix — floor(log2(size(x,1)))

  • vector — floor(log2(length(x)))

If numLevel is not specified, mlpt uses the maximum value.

Data Types: double

Name-Value Arguments

collapse all

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: PrimalMoments=3 computes a transform using three primal vanishing moments.

Number of dual vanishing moments in the lifting scheme, specified as one of the following: 2, 3 or 4.

Data Types: double

Number of primal vanishing moments in the lifting scheme, specified as one of the following: 2, 3, or 4.

Data Types: double

Prefilter before mlpt operation, specified as one of the following: "Haar" [1], "UnbalancedHaar", or "None".

Output Arguments

collapse all

MLPT coefficients, returned as a vector or matrix of coefficients, depending on the level to which the transform is calculated. coefs contains the approximation and detail coefficients.

Data Types: double

Sampling instants corresponding to output, returned as a vector or duration array of sample times obtained from x and t. The imlpt function requires T as an input. If the input t is a datetime or duration array, t is converted to units that allow for the stable computation of the mlpt and imlpt. Then T is returned as a duration array.

Data Types: double | duration

Coefficients per resolution level, returned as a vector containing the number of coefficients at each resolution level in coefs. The elements of coefsPerLevel are organized as follows:

  • coefsPerLevel(1) — Number of approximation coefficients at the coarsest resolution level.

  • coefsPerLevel(i) — Number of detail coefficients at resolution level i, where i = numLevel – i + 2 for i = 2, ..., numLevel + 1.

The smaller the index i, the lower the resolution. The MLPT is two times redundant in the number of detail coefficients, but not in the number of approximation coefficients.

Data Types: double

Scaling function moments, returned as a length(coefs)-by-P matrix, where P is the number of primal moments specified by the PrimalMoments name-value pair.

Data Types: double

Algorithms

Maarten Jansen developed the theoretical foundation of the multiscale local polynomial transform (MLPT) and algorithms for its efficient computation [1][2][3]. The MLPT uses a lifting scheme, wherein a kernel function smooths fine-scale coefficients with a given bandwidth to obtain the coarser resolution coefficients. The mlpt function uses only local polynomial interpolation, but the technique developed by Jansen is more general and admits many other kernel types with adjustable bandwidths [2].

References

[1] Jansen, Maarten. “Multiscale Local Polynomial Smoothing in a Lifted Pyramid for Non-Equispaced Data.” IEEE Transactions on Signal Processing 61, no. 3 (February 2013): 545–55. https://doi.org/10.1109/TSP.2012.2225059.

[2] Jansen, Maarten, and Mohamed Amghar. “Multiscale Local Polynomial Decompositions Using Bandwidths as Scales.” Statistics and Computing 27, no. 5 (September 2017): 1383–99. https://doi.org/10.1007/s11222-016-9692-8.

[3] Jansen, Maarten, and Patrick Oonincx. Second Generation Wavelets and Applications. London ; New York: Springer, 2005.

Version History

Introduced in R2017a

See Also

| |

Topics