MATLAB Wigner plot for Matching Pursuit atoms

Question

Using MATLAB I apply Matching Pursuit to approximate a signal. My problem is that I struggle to visualize the time-frequency representation of the selected atoms. I'm trying to produce a Wigner plot similar to the following image (source).

I have looked into the Wavelet Toolbox, Signal Processing Toolbox as well as the open source Time-Frequency Toolbox, but I'm possibly just using the wrong parameters, since my experience with signal processing is quite limited.

Example

Using this data my goal is to reproduce the plot from above.

% fit the signal using MP
itermax = 50;
signal = load('signal.txt');
dict = wmpdictionary(length(signal));
[signal_fit, r, coeff, iopt, qual, X] = wmpalg('OMP', signal, dict, ...
                                               'itermax', itermax);

% wigner plot of the simulated signal
tfrwv(signal_fit)  % wigner-ville function from time-frequency toolbox

% wigner plot of each atom
atoms = full(dict(:, iopt))  % selected atoms
for i = 1:itermax
    tfrwv(atoms(:, i))
end

Unfortunately, none of the resulting plots comes close to the target visualization. Note, that in the example I use tfrwv with standard parameters which I tweak with the GUI that it opens.

I'd greatly appreciate your help.

Update

I think I have now understood that one needs to use Gabor atoms to get blobs with shapes resembling stretched gaussians. Unfortunately, there are no Gabor functions in the predefined dicts of the Signal Processing Toolbox. However, this question helped me in implementing the needed dictionaries, such that I get atoms which look quite similar to the example:

Since my plots come close but are not perfect, there are still two questions open:

• Can all of the blobs that we see in the first example be modeled by Gabor atoms alone, or do I need another dictionary of functions?
• How can I combine the indidividual imagesc plots into a single visualization?

Answer 1

To answer your second question 'How can I combine the indidividual imagesc plots into a single visualization?'

If you have multiple 2d matrices that you want to superimpose and display using imagesc, I would suggest taking the element-wise maximum.

For example, I generate two 31x31 grids with gaussians with different mean and variance.

function F = generate2dGauss(mu, Sigma)
    x1 = -3:.2:3; x2 = -3:.2:3;
    [X1,X2] = meshgrid(x1,x2);
    F = mvnpdf([X1(:) X2(:)],mu,Sigma);
    F = reshape(F,length(x2),length(x1));
end

F1 = generate2dGauss([1 1], [.25 .3; .3 1]);
F2 = generate2dGauss([-1 -1], [.1 .1; .1 1]);

I can plot them with subplots as in your example,

figure; 
subplot(1,2,1);
title('Atom 1');
imagesc(F1);

subplot(1,2,2);
title('Atom 2');
imagesc(F2);

Or I can plot the per element maximum of the two grids.

figure;
title('Both Atoms');
imagesc(max(F1, F2));

You can also experiment with element-wise means, sums, etc, but based on the example you give, I think maximum will give you the cleanest looking result.

Possible pros and cons of different functions:

1. Maximum will work best if your atoms always have zero-valued backgrounds and no negative values. If the background is zero-valued, but the atoms also contain negative values, the negative values may be covered up by the background of other atoms. If your atom's overlap, the higher value will of course dominate.
2. Mean will make your peaks less high, but may be more intuitive where you have overlap between atoms.
3. Sum will make overlapping areas larger valued.
4. If you have non-zero backgrounds, you could also try using logical indexing. You would have to make some decisions about what to do in overlapping areas, but it would make it easy to filter out backgrounds.

Answer 2

Q. How can I combine the indidividual imagesc plots into a single visualization?

A. Use subplot to draw multiple plots, find below sample with 2 by 2 plots in a figure. Change your equations in code

x = linspace(-5,5);
y1 = sin(x);
subplot(2,2,1)
plot(x,y1)
title('First subplot')

y2 = sin(2*x);
subplot(2,2,2)
plot(x,y2)
title('Second subplot')

y3 = sin(4*x);
subplot(2,2,3)
plot(x,y3)
title('Third subplot')

y4 = sin(6*x);
subplot(2,2,4)
plot(x,y4)
title('Fourth subplot')