Strangely large value for 'abs(tf)' output from timefreq() function

[arnodelorme]

Makoto 2017-08-16 10:15:40 PDT Dear Andreas,

I would suggest to change units to "per sample" for time axis and frequency to normalized (implicitly).

This makes sense to me.

By the way, the title of this report says abs(tf) but if 'tf' is complex it should be conj(tf) to compute power... oops I made a mistake.

Makoto

[reply] [−]Comment 2Andreas Widmann 2017-08-16 04:08:51 PDT

This is still on my agenda. I will try to submit a patch in September.

The suggestion quoted below is a pragmatic solution, however, in two years nobody will understand why this was implemented this way. I would rather suggest to implement a more generic solution.

To my understanding the root of the problem is an incorrect mixture of units. The original Tallon-Baudry equations are based on normalized frequency units (as basically any text book equation). The current implementation now uses Hz and seconds for time axis but does not include sampling rate into the equation.

I would suggest to change units to "per sample" for time axis and frequency to normalized (implicitly). This way equations would be consistent with the literature. But just including sampling rate into the equation would be equally valid.

Best, Andreas

[reply] [−]Comment 1Arnaud Delorme 2017-08-07 11:48:40 PDT

Message from Andreas

Implement this in dftfilt3 fxArray_unit = fxArray * sqrt( pnts / max( sum( fxArray .^ 2 ) ) );

---

Dear all,

for me the discussion has not yet come to a satisfying conclusion. I would like to argue that it is essential that we agree on a common terminology. I expect that some (normalization) issues will automagically disappear (but as always some new issues come up ;)

(1) There are two families of transformations involved and intermixed in the discussion: Wavelet vs. STFT

In wavelet transforms time-frequency decomposition is performed with template waveforms (the mother wavelet) which are scaled (dilated or compressed) and translated. Due to scaling the time and frequency resolution changes across the time-frequency-plane (TF-plane): high frequency and low time resolution at low frequencies and low frequency and high time resolution at high frequencies. The wavelet transform is a multiresolution analysis (MRA) and this is not a side effect but the core feature. There exist various wavelets with different properties; Morlet wavelets are sine based and thus closely related to STFT and FT. However, by scaling the number of cycles in a Morlet wavelet is never changed. I have never read a definition of (Morlet) wavelet transform in a textbook changing cycles instead of scaling a wavelet. If you find one, please, send me the reference.

In contrast to the wavelet transform the STFT family has a constant fixed time and frequency resolution across the whole TF plane. There are various versions, in particular application of different window types but also variations of other parameters. As in Morlet wavelets a sinusoid is multiplied by a window function. To achieve constant frequency resolution of the atoms the number of cycles has to change with frequency (not arbitrarily but in a fixed relation as already annotated by Norman) while window width is kept constant.

There are time-frequency transformations not clearly belonging to any of these two categories as "multiresolution STFT" and others. Neither wavelet transform nor STFT is to be preferred in general. They are equally valid but targeting different applications. MRA is preferred if frequency resolution is considered more important at low frequencies and time resolution at high frequencies as e.g. in audition. Googling „wavelet vs. STFT“ gives dozens of presentations and papers comparing both methods including their (dis-)advantages and applications.

(2) Normalization has different requirements for wavelet transforms vs. STFT

Scaling wavelets changes their energy. This is disadvantageous as now the energy of the output of the wavelet transform would change with scale. Therefore, most commonly sets of (scaled) wavelets are normalized to have all the same energy. It not necessarily matters whether this is unit energy or not. But personally I would prefer unit energy for two reasons. First, because it gives a straight-forward and intuitively accessible interpretation of the output amplitudes: the amplitude of a sample in the output reflects the integrated energy of the input in the surface area covered by the wavelet in the TF plane at the given scale and delay weighted by the window function (also see my previous message how to compute the wavelet transform amplitude of a sinusoid). Second, because data and amplitudes can be more easily compared across studies and methods (also see example code below).

The atoms of the STFT, however, are already energy normalized "by design" as their window width is constant. Thus, energy normalization is not a core issue in STFT. Unit energy normalization might be advantageous for the reasons stated above but it will work also without. If one agrees to the above differentiation I would consider the fourth version in the code contributed by Mike ("Define empirical FWHM of wavelet in the frequency domain“) rather as a member of the STFT family (time and frequency resolution is constant). I have added a single set of these (once with and once without unit energy normalization) to my previous code hopefully demonstrating the differences between transforms and normalizations.

(3) Suggestions for EEGLAB

I would suggest to modify dftfilt3 so that by defaults wavelets are scaled to unit energy (as indirectly suggested by the documentation) and adding an option for backward compatible amplitude scaling. I would not add Gabor transform equivalent normalization to the multiresolution Morlet wavelet code (in case this is really needed this is a one-liner on the command line). Rather in my opinion a nice option would be to separately implement a textbook Gabor transform STFT in a separate function (actually reflecting Mike’s suggestion with minor changes; see Norman’s remarks).

(4) General remark

Sorry, I cannot withhold a more general „philosophical“ remark. Imho, the highest priority is replicability. Whatever analysis one performs should be described most thoroughly and detailed allowing for perfect replication and comparison of results across datasets, studies, researchers, labs, etc. This is considerably easier using textbook methods and terminology. For new methods and own approaches it is essential to describe them in a sufficiently detailed and accessible way. In my view this is the biggest problem with past and present literature applying TF transforms to electrophysiological data.

Best, Andreas

fs = 256; T = 4; pnts = T * fs; t = ( 0:pnts - 1 ) / fs; cycles = 7; F = 4:2:70;

% White noise rng( 0 ) signal = randn( fs * T, 1 ); % f = [ 10 50 ]; % signal = sin( 2 * pi * f( 1 ) * t )' + sin( 2 * pi * f( 2 ) * t )';

%% Morlet wavelet freqs = F * ( pnts / fs ); sigma_f = repmat( freqs / cycles, [ pnts 1 ] ); f = repmat( 0:pnts - 1, [ size( freqs, 2 ) 1 ] )' - repmat( freqs, [ pnts 1 ] ); cmorlf = exp( -f .^ 2 ./ ( 2 * sigma_f .^ 2 ) );

figure

% Analog to Gabor transform A = 2; cmorlf_gabor = A .* cmorlf;

tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf_gabor );

h = subplot( 4, 3, 3, 'YDir', 'normal' ); imagesc( t, F, abs( tf' ) .^ 2 ) colorbar set( h, 'YDir', 'normal' ); subplot( 4, 3, 2 ); plot( F, mean( abs( tf ) .^ 2 ) ) subplot( 4, 3, 1 ) plot( (0:pnts / 2) * fs / pnts, cmorlf_gabor( 1:pnts / 2 + 1, : ) )

% Unit energy A = sqrt( pnts ./ ( sqrt( pi ) * sigma_f ) ); cmorlf_energy = A .* cmorlf;

tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf_energy );

h = subplot( 4, 3, 6, 'YDir', 'normal' ); imagesc( t, F, abs( tf' ) .^ 2 ) colorbar set( h, 'YDir', 'normal' ); subplot( 4, 3, 5 ); plot( F, mean( abs( tf ) .^ 2 ) ) subplot( 4, 3, 4 ) plot( (0:pnts / 2) * fs / pnts, cmorlf_energy( 1:pnts / 2 + 1, : ) )

%% STFT % adapted from https://sccn.ucsd.edu/pipermail/eeglablist/2016/011536.html % by mikexcohen at gmail.com hz = (0:pnts - 1) * fs / pnts; fwhm = 10; for iF = 1:length( F ) s = fwhm*(2pi-1)/(4pi); % normalized width x = hz-F( iF ); % shifted frequencies fx = exp(-.5*(x/s).^2)'; % gaussian fxArray( :, iF ) = fx./max(fx); % gain-normalize end

tf = sqrt( 2 ) * ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* fxArray );

h = subplot( 4, 3, 9, 'YDir', 'normal' ); imagesc( t, F, abs( tf' ) .^ 2 ) colorbar set( h, 'YDir', 'normal' ); subplot( 4, 3, 8 ); plot( F, mean( abs( tf ) .^ 2 ) ) subplot( 4, 3, 7 ) plot( (0:pnts / 2) * fs / pnts, fxArray( 1:pnts / 2 + 1, : ) )

% Unit energy fxArray_unit = fxArray * sqrt( pnts / max( sum( fxArray .^ 2 ) ) );

tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* fxArray_unit );

h = subplot( 4, 3, 12, 'YDir', 'normal' ); imagesc( t, F, abs( tf' ) .^ 2 ) colorbar set( h, 'YDir', 'normal' ); subplot( 4, 3, 11 ); plot( F, mean( abs( tf ) .^ 2 ) ) subplot( 4, 3, 10 ) plot( (0:pnts / 2) * fs / pnts, fxArray_unit( 1:pnts / 2 + 1, : ) )

[reply] [−]DescriptionMakoto 2016-08-08 18:09:57 PDT

Hi this is Makoto. This issue has been there since timefreq() was written. To sum, this is because sum((wavelet) x (data)) is not normalized by length of wavelet... that's my intuition. This is why abs(tf) is extremely large and changes a lot when wavelet cycles are changed. See below for original communication on EEGLAB mailing list with Niko Busch.

Makoto

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% On Wed, Aug 3, 2016 at 10:25 AM, Niko Busch [email protected] wrote: Hi everyone,

can you help me understand the output of the timefreq function? The main output is "tf", and abs(tf) is a measure of spectral power. But what exactly do the units of tf represent?

Below, I attach a simple script that simulates a synthetic oscillation of known frequency and amplitude. Using timefreq on this signal, I would expect the values of abs(tf) to match the amplitude of that oscillation. But as you see, the values are orders of magnitude larger than the amplitude of the oscillation. Furthermore, I notice that the magnitude of the values in tf depends on the sampling frequency of the input signal and the number of wavelet cycles.

So, what do values in tf represent and how should I "normalize" these values to match the amplitude of the input signal?

Thanks, Niko

%% Create sine wave clear all D = 4; % total signal duration in seconds. sigD = 1; % duration of the test oscillation within the signal. F = 10; % frequency of the test oscillationin Hz. P = .25; % Phase of the test oscillation. 2 pi radians = 360 degrees srate = 256; % sampling rate, i.e. N points per sec used to represent sine wave. T = 1/srate; % sampling period, i.e. for this e.g. points at 1 ms intervals t = [T:T:D]; % time vector.

sigpoints = length(t)/2 - (sigDsrate)/2:(length(t)/2 + (sigDsrate)/2)-1; mysig = zeros(1,Dsrate); mysig(sigpoints) = sin(2Ft(sigpoints)pi+ 2piP);

%% TF computation [tf, outfreqs, outtimes] = timefreq(mysig', srate, ... 'cycles', 4, 'wletmethod', 'dftfilt3', 'freqscale', 'linear', ... 'freqs', F);

%% Plot figure; hold all plot(t,mysig); plot(outtimes./1000,abs(tf)) xlabel('Time (seconds)'); ylabel('Amplitude'); legend('input signal', 'wavelet result')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Niko, I know your question is valid because I have been wondering this issue over a decade! Actually, the more critical question is why the raw data before dB conversion is so large (such as several thousands, when data are only +/- 1uV peak-to-peak). It's time for me to find it out myself.

I investigated it running your example code. I found that in timefreq line 453-454,

tmpX = sum(wav .* tmpX); tmpall( freqind, index, :) = tmpX;

'wav' is a wavelet kernel generated by dftfilt3 (default) and tmpX is the chunk of EEG data provided, and the output tmpX is a complex number, and this is eventually passed to 'tf' which is the final complex output.

When I increased 'cycles' from 4 to 6 and 8, the length of 'wav' were 115, 171, 229 (by the way this can be calculated by 'wavelet = dftfilt3(freqs, [4|6|8], srate); length(wavelet)), the value of 'tmpX' increased as well--this is so obvious, because 'tmpX' is a simple sum of wav.*tmpX without normalization, so as 'wav' becomes longer, sum of 'tmpX' becomes larger.

I thought that the output 'tmpX' needs to be normalized by the length of data, such as /length(wav) to allow users to compare raw output. So far this was not a problem because most of us have been using dB conversion so scale was not a problem. However, as Niko is asking, it is confusing that the non-normalized 'tf' behaves in non-intuitively.

Makoto