Wavelet Multiscale Principal Components Analysis

This section demonstrates the features of multiscale principal components analysis provided in the Wavelet Toolbox™ software.

The aim of multiscale PCA is to reconstruct, starting from a multivariate signal and using a simple representation at each resolution level, a simplified multivariate signal. The multiscale principal components generalizes the normal PCA of a multivariate signal represented as a matrix by performing a PCA on the matrices of details of different levels simultaneously. A PCA is also performed on the coarser approximation coefficients matrix in the wavelet domain as well as on the final reconstructed matrix. By selecting the numbers of retained principal components, interesting simplified signals can be reconstructed.

This example uses noisy test signals. In this section, you will:

Load a multivariate signal.
Perform a simple multiscale PCA.
Display the original and simplified signals.
Improve the obtained result by retaining less principal components.

Load a multivariate signal by typing at the MATLAB^® prompt:
```
load ex4mwden  
whos
```
Name Size Bytes Class
covar 4x4 128 double array
x 1024x4 32768 double array
x_orig 1024x4 32768 double array
The data stored in matrix x comes from two test signals, Blocks and HeavySine, and from their sum and difference, to which multivariate Gaussian white noise has been added.
Perform a simple multiscale PCA.
The multiscale PCA combines noncentered PCA on approximations and details in the wavelet domain and a final PCA. At each level, the most significant principal components are selected.
First, set the wavelet parameters:
```
level= 5; 
wname = 'sym4';
```
Then, automatically select the number of retained principal components using Kaiser's rule by typing
```
npc = 'kais';
```
Finally, perform multiscale PCA:
```
[x_sim, qual, npc] = wmspca(x ,level, wname, npc); 
```

Name	Size	Bytes	Class
`covar`	`4x4`	`128`	`double array`
`x`	`1024x4`	`32768`	`double array`
`x_orig`	`1024x4`	`32768`	`double array`

Display the original and simplified signals:

kp = 0; 
for i = 1:4  
    subplot(4,2,kp+1), plot(x (:,i)); set(gca,'xtick',[]);
    axis tight; 
    title(['Original signal ',num2str(i)]) 
    subplot(4,2,kp+2), plot(x_sim(:,i)); set(gca,'xtick',[]);
    axis tight; 
    title(['Simplified signal ',num2str(i)]) 
    kp = kp + 2;
end

The results from a compression perspective are good. The percentages reflecting the quality of column reconstructions given by the relative mean square errors are close to 100%.

qual

qual =  

   98.0545   93.2807   97.1172   98.8603

Improve the first result by retaining fewer principal components.
The results can be improved by suppressing noise, because the details at levels 1 to 3 are composed essentially of noise with small contributions from the signal. Removing the noise leads to a crude, but large, denoising effect.
The output argument npc contains the numbers of retained principal components selected by Kaiser's rule:
```
npc  

npc =   
     1     1     1     1     1     2     2
```
For d from 1 to 5, npc(d) is the number of retained noncentered principal components (PCs) for details at level d. The number of retained noncentered PCs for approximations at level 5 is npc(6), and npc(7) is the number of retained PCs for final PCA after wavelet reconstruction. As expected, the rule keeps two principal components, both for the PCA approximations and the final PCA, but one principal component is kept for details at each level.
To suppress the details at levels 1 to 3, update the npc argument as follows:
```
npc(1:3) = zeros(1,3);

npc  

npc = 
0     0     0     1     1     2     2
```
Then, perform multiscale PCA again:
```
[x_sim, qual, npc] = wmspca(x, level, wname, npc); 
```

Display the original and final simplified signals:

kp = 0; 
for i = 1:4  
    subplot(4,2,kp+1), plot(x (:,i)); set(gca,'xtick',[]);
    axis tight; 
    title(['Original signal ',num2str(i)]); set(gca,'xtick',[]);
    axis tight;  
    subplot(4,2,kp+2), plot(x_sim(:,i)); set(gca,'xtick',[]);
    axis tight;  
    title(['Simplified signal ',num2str(i)]) 
    kp = kp + 2;
end

As shown, the results are improved.