We have decided to get to the bottom of the unusual behavior of MPTK,
since the next steps of our work on CMP depend on it.
This entails comparing the results from my MATLAB implementation with those of MPTK (and CMPTK) on the same dictionary.
I have decomposed one signal of dimension 1024 samples,
with a dictionary of modulated Gaussian windows of scale 64 samples,
with a hop size of 8 samples.
The MPTK dictionary is defined as follows (couldn't use pre html tags for some reason, so I attach it as a png):
This "windowOpt = 0" business means that the variance of the Gaussian window is set to the default 0.02.
In my MATLAB code I have attempted to do the same. Here is how MPTK computes the window (actually only half of it is needed) of even length length:
optional = 0.02;
optional = 1/(2*optional*(length+1)*(length+1));
for ( i = 0; i < length/2; i++, p1++, p2-- ) {
k = (double) i - ( (double) (length - 1) )/2.0;
newPoint = exp(-k*k*optional);
*p2 = *p1 = (Dsp_Win_t) newPoint;
}
Put together, p1 and p2 define the entire window.
Here is how I compute the same Gaussian window in MATLAB:
winVals = zeros(1,size);
optional = 0.02;
optional = 1/(2*optional*(size+1)*(size+1));
counter = 1;
for ii=0:floor(size/2)-1
kk = ii - (size-1)/2;
winVals(counter) = exp(-kk.^2*optional);
winVals(size - counter + 1) = exp(-kk.^2*optional);
counter = counter + 1;
end
Both dictionaries have 3993 complex atoms.
So far, everything checks out.
In the graph below, I show the ratios of residual energies due to my implementation to those produced by MPTK/CMPTK.
Above the periwinkle dashed line, the residual is larger in MATLAB that in MPTK; and vice versa below the line.
We see for the most part the differences of MP/MPTK are negligible.
This could be due to some indexes being one off, some of which I did find in my code and fixed.
For CMP/CMPTK however, the story is different. (I am only performing one cycle of refinement in each case.)
There are no reported increases of the residual energy in CMP,
but in CMPTK there are plenty, beginning with atom 11.
Let's have a look at the residual decreases of CMP in MATLAB:
2: Energy reduction of 0.27118%
3: Energy reduction of 0.13418%
4: Energy reduction of 0.37532%
5: Energy reduction of 0.073902%
6: Energy reduction of 15.5745%
7: Energy reduction of 5.8543%
8: Energy reduction of 1.0909%
9: Energy reduction of 1.764%
10: Energy reduction of 1.497%
And now of CMPTK
2: Energy reduction of 0.266262% from 1 of 1 possible cycles
3: Energy reduction of 0.131654% from 1 of 1 possible cycles
4: Energy reduction of 0.40187% from 1 of 1 possible cycles
5: Energy reduction of 0.0768828% from 1 of 1 possible cycles
6: Energy reduction of 4.58075% from 1 of 1 possible cycles <-----
7: Energy reduction of 11.134% from 1 of 1 possible cycles
8: Energy reduction of 2.4887% from 1 of 1 possible cycles
9: Energy reduction of 1.68054% from 1 of 1 possible cycles
10: Energy reduction of 1.72243% from 1 of 1 possible cycles
All looks compatible until iteration 5 (arrow). We see that in the graph above too, when CMPTK performance takes a nose dive like my recent investments in CMP, Inc.
So what is the atom that causes trouble at iteration 11?
There is nothing unusual about it. It begins at sample 80 (shouldn't that be 79 if we begin indexing at 0?), has a frequency of 0.015625*8192 = 128 Hz, an amp of 0.277851, and a phase of -1.82707. This atom causes the residual energy to increase 1.146%.
Here is where we take a look at the residual just before and just after this atom comes into play.
(Editors note: several hours passed between these two paragraphs, during which I had a long think to the bus stop, and some more think over dinner.)
Upon returning to home, I decided to perform the ultimate test of sanity: put two atoms from the dictionary in a signal, and see what happens. With this I believe I have cornered the problem. Take a look at the following movie.
My MATLAB implementation vs. MPTK. One is slow, and one is speedy.
Which one will win?
ex1_MP_atoms.mov. Clearly, MPTK has made a very poor choice. And I believe it all comes down to the position index selecting the atom in the block.
What happens when I make off by one in my MATLAB MP those positions at which I am searching for the largest magnitude projection (i.e., the atom hop size)?:ex1_MP_atoms2.mov.
Not only do I pick up the same wrong atoms as MPTK, my computer's battery nears empty. The 2-norm of the differences of the residuals is 1e-4 --- which is not too large, but not too small.
Finally, when I shift the dictionary such that the atoms line up with how MPTK expects them, and keep my erroneous MP, I get the following: ex1_MP_atoms3.mov. Yay! It works!!!
But wait. That's not all. Often I see the following behavior: ex1_MP_atoms4.mov. That first projection was a doozy for MPTK. It left a massive dirac. And those 3 atoms aren't even close to each other to interfere. And here is one where it happens twice: ex1_MP_atoms5.mov. My MP converges after 3 steps with an SRR of 94 dB; while MPTK has a paltry 31 dB SRR. The only thing I can think of is that MPTK is so fast that it becomes extremely excited when something fits so perfectly that it twists its atoms in knots before subtraction.
So, not only do we need to change the positions at which MPTK searches for maximums,
but also something in the atom creation causing those diracs. And then CMPTK will hopefully be ready for the next step.]]>
Now, with the corrections, I observe the decays. The "MPold" decay is that produced by the uncorrected MPTK.
"MP" shows that of the new code.
Only in Attack and Sine do we see much difference;
and at times in Sine the previous version of MPTK beats the corrected version. (Such is the behavior of greedy algorithms. I will write a Po'D about this soon.)
Anyhow, the decays of CMP-\(\ell\) (where the number denotes the largest number of possible cycles of refinement, but I suspend refinement cycles when energyAfter/energyBefore > 0.999), comports with the decays I see in my MATLAB implementation (see above).
So, now I am comfortable moving on.
Below we see the decays and cycle refinements for three different CMPs for these four signals. (Note the change in the y axes.) Bimodal appears to benefit the most in the short term from the refinement cycles, after which improvement is sporadic.
The modeling of Sine has a flurry of improvements.
It is interesting to note that as \(\ell\) increases, we do not necessarily see better models with respect to the residual energy. For instance, for Attack, the residual energy for CMP-1 beats the others.
And briefly back to the glockenspiel signal, below we see the decays and improvements
using a multiscale Gabor dictionary (up to atoms with scale 512 samples).
]]>
Below, I show the time-domain residual signals resulting from both MP and CMP-1.
For the most part, the CMP-1 error signal is below that of MP; but strangely the first attack causes more problems for CMP-1 than the other attacks.
