Forum: help

RE: Normalization of gabor atoms [ Reply ]
By: Supratim Ray on 2008-07-21 13:59

Dear Remi

To reduce the problem of a sudden discontinuity in the atoms (which are non-zero for windowLen followed by NFFT-windowLen zeros), I made the windowLength equal to the FFTsize, and decreased the value of the window-opt parameter. But now there appears to be a different problem: the dictionary is no longer complete.

In your implementation, the window is given by exp(-1/(2*opt*(N+1)^2)), where N is the window size. When compared to the standard form of exp(-pi/s^2), yields a scale of s such that s^2 = 2*pi*opt*(N+1)^2. Hence if opt is less than 1/(8*pi) ~0.04, s is less than N/2. The scale of a gaussian atom determines its spread (temporal spread of s and spectral spread of 1/s). It appears that for small opt values, the edge of the signal is not covered.

Consider a window and FFT size of N=512, with a Windowshift = 64. With only these atoms, the decomposition is just like doing a STFT with a constant window. The problem is that the maximum center position of an atom in your implementation is only upto L - N/2, where L is the length of the signal, which leaves the last N/2 points "uncovered" if the scale is less than N/2.

As an example, I tried to reconstruct a signal (that I sent in the previous email) with the following dictionary

<?xml version="1.0" encoding="ISO-8859-1"?>
<dict>
<libVersion>0.2</libVersion>
<blockproperties name="GAUSS-WINDOW">
<param name="windowtype" value="gauss"/>
<param name="windowopt" value="0.02"/>
</blockproperties>

<block uses="GAUSS-WINDOW">
<param name="type" value="gabor"/>
<param name="windowLen" value="512"/>
<param name="windowShift" value="64"/>
<param name="fftSize" value="512"/>
</block>
</dict>

The decomposition, not surprisingly, fails to pick out the atom centered at n=1000 completely. However, as you mentioned before, this is a problem only at the edges of the signal (not a big issue for very long signals I suppose). For short duration signals, however, this would lead to some problems.

regards,
Supratim

RE: Normalization of gabor atoms [ Reply ]
By: Supratim Ray on 2008-07-18 23:49

[forum:98166]

Dear Remi,

Thank you for your detailed reply. I spent a lot of time trying to understand the algorithm in detail, but I still have some questions. Please let me know what you think of these issues:

1. Normalization issue: As you said, the 'amp' parameter is not the inner product, since the atoms are not normalized. However, in this case the inner product is just the product of the amp value and the actual norm of the atom (which can be obtained from the reconstructed atoms). I plotted the product of amp and norm expecting to see a monotonic decrease with atom number, but in fact that does not happen. For atoms with very small temporal support or very low frequency, there are deviations.

Just to clarify, you are indeed doing the MP on real atoms (equation 55 of Mallat and Zhang), right? Not a complex MP (equation 42 of MZ) followed by some ways to make the atoms real? I'm sorry I looked at your compute_energy code, but could not follow some details. I'm sending a separate email with my data/results and some figures, since I can't attach them here. Kindly look at my attachment and let me know what you think.

2. Periodization issues: I did not follow your reasons for not periodizing the atoms. As far as I could understand, periodization is necessary because the FFT analysis assumes that the signals are periodic (the multiplication by a complex phase, for example, results in a cyclic shift in the FFT). In Mallat and Zhang, periodization was done to get an analytic solution of the inner products between atoms (Lemma 6, equation 121-122). If this is true, the "boundaries" I was talking about are not the actual end points of the signal - they are the end points between which FFT is taken. For the original M&Z MP, the atoms were always equal to the size of the signal, so the boundary issues were limited to the actual boundaries of the signal. I think the issue is more serious here because you are taking FFT of different sizes repeatedly. Further, the atoms on which the signal is projected are often discontinuous - they have a length equal to 'len' followed by (K - len) zeros, where K is the FFT size. Unless the window_opt parameter is very small, this will lead to discontinuities in the residue which may introduce broad-band noise.

I'm sorry I have so many questions. My main research interest is in studying brain signals and trying to understand the neural correlates of the so called "gamma" oscillations, which are highly non-stationary low amplitude signals thought to be correlated with higher brain functions. These are best picked out by short atoms with small temporal support, so I'm especially worried about their representation.

Thanks
Supratim

NB: I'm sending a doc file and some attachments separately.

RE: Normalization of gabor atoms [ Reply ]
By: Rémi Gribonval on 2008-07-16 08:02

[forum:96958]

Hi Supratim,

Thanks for your interesting and relevant question.
>First, is the 'amp' parameter in the atom
>structure (obtained from the function bookread)
>the value of the inner product between the atom
>and the signal?

No, it is not necessarily, although it always matches with the inner product for all currently existing atom types except Gabor and Harmonic.
However, at each stage of MP, the Gabor atom which is selected and removed is exactly the one with maximum correlation with the residual, as explained further below.

You should think of an atom as a description of a signal. The field 'amp' describes a multiplying factor applied to a waveform described by the other parameters, but the waveform is not necessarily normalized. To understand the meaning of 'amp' for Gabor atom you should refer to the help of MP_Gabor_Atom_Plugin_c::build_waveform, which is defined in gabor_atom_plugin.cpp and can also be found in the doxygen documentation.

As a matter of fact, when MP is run, the selection of a Gabor atom is done as follows
-the signal is multiplied by the unit energy window and FFT
-the complex FFT values correspond to inner products between the signal and complex Gabor atoms. These complex atoms are simply the product of the window by a complex exponential, and they are of unit energy.
-a real Gabor atom is a linear combination of two conjugate complex Gabor atoms (at +/- frequency), and therefore belongs to the linear span of these two atoms.
-the best real Gabor atom at a given frequency is the projection of the signal onto the linear span of the two conjugate complex atoms. The projection is computed using the complex FFT values at +/- frequency (which are conjugate to one another) as well as the cross correlation between the conjugate complex atoms. This is done in MP_Gabor_Block_c::compute_energy() where you will find additional documentation on this.

>Are there any theoretical issues we need to
>watch out because of limiting our basis signals
>to a finite range?
Not that I know of. MP can decompose a signal on an arbitrary dictionary of atoms. Choosing periodic atoms or finite support atoms amounts to choosing how the signal is modelled. Since our aim in MPTK was to allow the decomposition of very large signals, we did not provide code for dictionaries of periodized atoms supported on the whole signal support. If suppose periodized atoms could be added by writing the appropriate plugin, if this is needed.

I hope this helps!

Remi.

Normalization of gabor atoms [ Reply ]
By: Supratim Ray on 2008-07-15 20:23

[forum:96657]

Hi,

Sorry for spamming here, I'm running into a lot of problems while trying to understand the MPTK code.

First, is the 'amp' parameter in the atom structure (obtained from the function bookread) the value of the inner product between the atom and the signal? I wrote a reconstruction program in matlab to generate the gabor basis signals from the paramters in the 'atom' structure (code attached at the bottom of this message), but found that the reconstruction works only if we do not normalize the atoms. If we multiply these basis signals with the 'amp' parameter, we get the reconstructed signal perfectly (I compared the matlab reconstruction with the reconstruction from the mpr command), but without normalization of the gabor basis signals the amp parameter is not equal to the inner product.

Looking back into your code (dsp_windows.c), I see that only the window function is normalized (line 398 to 407). This works well for signals with long temporal support (large values of atom.len), for which the entire signal has a norm of about 0.5. For short transient signals, the norms begin to deviate from 0.5, which is understandable, because for short signals the shape of the signal is dependent on the phase of the cosine part, and does not always integrate to 0.5.

This point is important because the algorithm chooses atoms based on the value of their inner product, and this might bias the selection of atoms with short temporal support.

Second, I see that the reconstruction works only when the signal is considered zero outside the 'len' parameter, which is consistent with your concept of atoms as short signals with an existence of their own. However, if the window_opt_parameter is not chosen to be very small, this introduces 'sharp edges' in the signal after subtraction of an atom. This is in contrast to the original MZ code where the basis functions were periodized over N points (equation 63, similar to the 'wrapping' option I mentioned in a previous message). Are there any theoretical issues we need to watch out because of limiting our basis signals to a finite range?

Sorry for so many questions.
Thank you
Supi

% Matlab code for reconstruction of gabor atoms

% This program reconstructs the atom given its description (the structure
% 'atomDesc'), signalLength and sampleRate (Fs). The strcuture atom is obtained from
% the book, and follows the conventions used in MPTK.

function [atom] = reconstructAtom(atomDesc,signalLength,Fs)

twopi = 6.28318530717958647692;

if ~exist('signalLength') signalLength = 1024; end
if ~exist('Fs') Fs=1000; end

% Properties of the atom
N = atomDesc.len; % temporal support of the atom
frontZeros = zeros(1,atomDesc.pos);
backZeros = zeros(1,signalLength - (atomDesc.pos+atomDesc.len));

if strcmp(atomDesc.type,'gabor') % Gabor

% Generate the window, as done in dsp_windows.c program
opt = 1/(2*atomDesc.windowOpt*(N+1)*(N+1));
winArray = (0:N-1) - (N-1)/2;
window = exp(-opt*winArray.^2);
normalizedWindow = window/sqrt(sum(window.*window)); % normalize window
cosine_part = cos(twopi*atomDesc.freq*(0:N-1)+atomDesc.phase);

atom = cosine_part.*normalizedWindow;
end

atom = [frontZeros atom backZeros];