Multiscale principal component analysis
returns a simplified version
pca_params] = wmspca(
xsim of the input matrix
x obtained from the wavelet-based multiscale principal component analysis
(PCA). The wavelet decomposition is performed using the decomposition level
level and the wavelet
Wavelet Principal Component Analysis of Noisy Multivariate Signal
Use wavelet multiscale principal component analysis to denoise a multivariate signal.
Load the dataset consisting of four signals of length 1024. Plot the original signals and the signals with additive noise.
load ex4mwden for i = 0:3 subplot(4,2,2*i+1) plot(x_orig(:,i+1)) axis tight title(['Original signal ',num2str(i+1)]) subplot(4,2,2*i+2) plot(x(:,i+1)) axis tight title(['Noisy signal ',num2str(i+1)]) end
Perform the first multiscale wavelet PCA using the Daubechies least-asymmetric wavelet with four vanishing moments,
sym4. Obtain the multiresolution decomposition down to level 5. Use the heuristic rule to decide how many principal components to retain.
level = 5; wname = 'sym4'; npc = 'heur'; [x_sim,qual,npcA] = wmspca(x,level,wname,npc);
Plot the result and examine the quality of the approximation.
for i = 0:3 subplot(4,2,2*i+1) plot(x(:,i+1)) axis tight title(['Noisy signal ',num2str(i+1)]) subplot(4,2,2*i+2) plot(x_sim(:,i+1)) axis tight title(['First PCA ',num2str(i+1)]) end
qual = 1×4 97.4372 94.5520 97.7362 99.5219
The quality results are all close to 100%. The
npc vector gives the number of principal components retained at each level.
Suppress the noise by removing the principal components at levels 1-3. Perform the multiscale PCA again.
npcA(1:3) = zeros(1,3); [x_sim,qual,npcB] = wmspca(x,level,wname,npcA);
Plot the result.
for i = 0:3 subplot(4,2,2*i+1) plot(x(:,i+1)) axis tight title(['Noisy signal ',num2str(i+1)]) subplot(4,2,2*i+2) plot(x_sim(:,i+1)) axis tight title(['Second PCA ',num2str(i+1)]) end
x — Multisignal
Multisignal, specified as a real-valued matrix. The matrix
contains P signals of length N stored column-wise
(N > P).
level — Level of decomposition
Level of decomposition, specified as a positive integer.
not enforce a maximum level restriction. Use
wmaxlev to ensure that the wavelet coefficients are free from boundary effects.
If boundary effects are not a concern, a good rule is to set
than or equal to
N is the signal length.
wname — Analyzing wavelet
character vector | string scalar
Analyzing wavelet, specified as a character vector or string scalar. The wavelet must be
orthogonal or biorthogonal. Orthogonal and biorthogonal wavelets are designated as type 1 and
type 2 wavelets respectively in the wavelet manager,
Valid built-in orthogonal wavelet families begin with
'symN', where N is the number of vanishing moments for all families except
fk, N is the number of filter coefficients.
Valid biorthogonal wavelet families begin with
'rbioNd.Nr', where Nr and Nd are the number of vanishing moments in the reconstruction (synthesis) and decomposition (analysis) wavelet.
Determine valid values for the vanishing moments by using
waveinfo with the wavelet family short name. For example, enter
wavemngr('type',WNAME) to determine if a wavelet is orthogonal (returns 1)
or biorthogonal (returns 2).
npc_in — Principal components parameter
Principal components parameter, specified as a vector, character vector, or string scalar.
npc_inis a vector, then it must be of length
level+2. The vector
npc_incontains the number of retained principal components for each PCA performed:
npc_in(d)is the number of retained noncentered principal components for details at level
d, for 1 ≤
npc_in(level+1)is the number of retained non-centered principal components for approximations at level
npc_in(level+2)is the number of retained principal components for final PCA after wavelet reconstruction.
npc_inmust be such that 0 ≤
npc_in(d)≤ P, where P is the number of signals, for 1 ≤
"kais", then the number of retained principal components is selected automatically using Kaiser's rule. Kaiser's rule keeps the components associated with eigenvalues exceeding the mean of all eigenvalues.
"heur", then the number of retained principal components is selected automatically using the heuristic rule. The heuristic rule keeps the components associated with eigenvalues greater than 0.05 times the sum of all eigenvalues.
"nodet", then the details are "killed" and all the approximations are retained.
extmode — Extension mode
'spd' | ...
Extension mode used when performing the wavelet decomposition, specified as:
DWT Extension Mode
Smooth extension of order 0
Smooth extension of order 1
Symmetric extension (half point): boundary value symmetric replication
Symmetric extension (whole point): boundary value symmetric replication
Antisymmetric extension (half point): boundary value antisymmetric replication
Antisymmetric extension (whole point): boundary value antisymmetric replication
If the signal length is odd and
dec — Wavelet decomposition structure
xsim — Simplified multivariate multisignal
Simplified multivariate multisignal, returned as a matrix. The dimensions of
xsim equal the dimensions of
qual — Quality of column reconstructions
Quality of column reconstructions, returned as a vector of length P,
where P is equal to
contains the quality of column reconstructions given by the relative mean square errors in
npc_out — Number of retained principal components
Number of retained principal components, returned as a vector. If
npc_in is a vector, then
decsim — Wavelet decomposition
Wavelet decomposition of the simplified multisignal
as a structure with the following fields:
'c'(column), indicator of decomposition direction
level— Level of wavelet decomposition
wname— Wavelet name
dwtFilters— Structure with four fields:
LoD— Lowpass decomposition filter
HiD— Highpass decomposition filter
LoR— Lowpass reconstruction filter
HiR— Highpass reconstruction filter
dwtEXTM— DWT extension mode
dwtShift— DWT shift parameter (0 or 1)
dataSize— Size of
ca— Approximation coefficients at level
cd— Cell array of detail coefficients, from level 1 to level
k from 1 to
level, are matrices, where the
coefficients are stored as columns.
pca_params — PCA parameters
PCA parameters, returned as a structure array of length
pca_params(d).pcis a P-by-P matrix of principal components. The columns are stored in descending order of the variances.
pca_params(d).variancesis the principal component variances vector.
pca_params(d).npc = npc_out
The multiscale principal components generalizes the usual PCA of a multivariate signal seen as a matrix by performing simultaneously a PCA on the matrices of details of different levels. In addition, a PCA is performed also on the coarser approximation coefficients matrix in the wavelet domain as well as on the final reconstructed matrix. By selecting conveniently the numbers of retained principal components, interesting simplified signals can be reconstructed.
 Bakshi, Bhavik R. “Multiscale PCA with Application to Multivariate Statistical Process Monitoring.” AIChE Journal 44, no. 7 (July 1998): 1596–1610. https://doi.org/10.1002/aic.690440712.