newtimef() does not scale FFT or Wavelet Power Spectra correctly

[arnodelorme]

This an original report on Bugzilla

When newtimef.m computes power spectra using either FFT or Morlet Wavelets, the results are scaled incorrectly. This has no effect on computed ERSP or ITC results when baseline correction is applied. However, it affects the returned value of 'powbase' (the mean power in the baseline interval), and if baseline correction is turned off, it affects the returned value of 'ersp'. The reported power changes inappropriately with parameters such as sample rate and window size.

Related BUG reports: #446, #1032, #661, #1920

EEGLAB version: 14_1_1b

Details:

if newtimef() is called with 'baseline'=NaN then baseline scaling is disabled, and the function produces a time-freq plot showing power in units of either uV^2/Hz or dB (relative to 1uV^2/Hz) depending upon the 'scale' parameter. These units imply that the function should be computing 'power spectral density' which, in my opinion, is the correct option when examining broad-band signals such as EEG data.

Critically, estimates of total power in the time-domain version of a signal should equal estimates of total power in the frequency-domain representation of that signal. i.e. the variance of the time-domain signal should equal the integral (area under the curve) of the power spectral density. This relationship holds theoretically for FFT analysis, and it would be conventional to scale the results from wavelet analysis in a comparable fashion. For example: for a sine wave with peak amplitude A, the total power (variance) of the time-domain signal should be A*A/2, irrespective of sample rate or analysis procedure.

Furthermore, in addition to the power spectral distribution (P), newtimeF() also returns 'alltfX' which contains a complex array PROPORTIONAL to the one-sided amplitude spectrum. At numerous locations within EEGLAB, power estimates are derived from alltfX using the relationship "P = alltfX .* conj(alltfX)". alltfX should be scaled so that this holds for both FFT and wavelet analysis.

Separate changes are required for the FFT analysis and the wavelet analysis:

FOR THE FFT ANALYSIS:

Within newtimef, the FFT calculations are performed by calling timefreq(). After this, the FFT results are scaled and power values computed. The relevant lines of code are 1185-1190:

alltfX = 2/0.375*alltfX/g.winsize;   % TF and MC (12/11/2006): normalization, divide by g.winsize
P  = alltfX.*conj(alltfX);           % power  

% TF and MC (12/14/2006): multiply by 2 account for negative frequencies,
% and ounteract the reduction by a factor 0.375 that occurs as a result of 
% cosine (Hann) tapering. Refer to Bug 446
% Modified again 04/29/2011 due to comment in bug 1032

I recommend changing these lines of code to:

FreqBinWidth = g.srate / (g.winsize * g.padratio);

alltfX = sqrt(2)/sqrt(FreqBinWidth)*sqrt(8/3)*sqrt(g.padratio)/(g.winsize*g.padratio)*alltfX;

P  = real( alltfX .* conj(alltfX) );

% Notes: alltfX returned by timefreq() is a complex array containing the
% lower half of the unnormalized two-sided amplitude spectrum from
% the FFT (for positive frequencies > 0 and <= Nyquist frequency).
%
% FreqBinWidth (in Hz) is the interval between frequencies extracted by 
% the FFT and is equal to the inverse of the PADDED window size (in
% seconds).
%
% (1) Normalise the amplitude spectrum in alltfX by dividing by the TOTAL
%     number of sample points in the zero-padded 
%     window (g.winsize * g.padratio).
% (2) Multiply alltfX by sqrt(g.padratio) to account for loss of power 
%     when zero padding.
% (3) Multiply alltfX by sqrt(8/3) to account for loss of power when 
%     applying hanning window.
% (4) Divide alltfX by sqrt(FreqBinWidth) so that computed power will 
%     be 'Power Spectral Density' in units of uV^2/Hz
% (5) Multiply alltfX by sqrt(2) to convert from two-sided to
%     one-sided spectrum.
%
% Note: alltfX now is a complex valued array that is PROPORTIONAL to 
% the discrete 1-sided amplitude spectrum, where:
%
% AmpSpectrum = sqrt(2) * sqrt(FreqBinWidth) * alltfX;
%
% AmpSpectrum contains the 1-sided discrete amplitude spectrum  
% representing amplitude and phase at each frequency extracted by the FFT,
% scaled to preserve the TOTAL POWER in the original signal.
% For this amplitude spectrum, each corresponding sine wave has a 
% peak amplitude (in uV) = abs(AmpSpectrum),
% and phase (in degrees) = angle(AmpSpectrum)/pi*180.
%
% (6) Convert alltfX to one-sided 'Power Spectral Density', P. 
%     Note that P is always a real quantity, so can take the real part
%     of P with no loss in inforation.
%
% P is a real array containing the Power Spectral Density of
% the signal in units uV^2/Hz. The INTEGRAL (area under the curve) of P
% equals the total power (variance) of the original signal.

ISSUES WITH THE ORIGINAL FFT CODE:

1. Adjusting for the hanning window involves multiplying the POWER spectrum by 1/.375 = 8/3. The original code incorrectly applies this scaling to the AMPLITUDE spectrum. The new code scales the amplitude spectrum by sqrt(8/3) so that both the amplitude and power spectrum are correct.

2. EEGLAB uses the formula "Power = Amplitude . conj(Amplitude)" to compute the power spectrum. This formula is correct for use with a TWO-SIDED amplitude spectrum with units of uV, or with a ONE-SIDED amplitude spectrum using units of uVrms (where 1 uV = 1/sqrt(2) uVrms). However, the original code computes the one-sided amplitude spectrum by multiplying the two-sided spectrum by 2. This produces a ONE-SIDED amplitude spectrum with units uV, for which the appropriate power formula would be: "Power = Amplitude . conj(Amplitude)/2".

3. The original code did not scale 'alltfX' correctly for padratio (even though the Power was unaffected).

4. The original code did not convert from 'power spectrum' to 'power spectral density'. Thus P had units uV^2, whereas it was being displayed when plotted with units uV^2/Hz. It needs to be power spectral density for consistency with the wavelet analysis.

TEST PROCEDURE:

Create a test signal consisting of a sine wave with peak amplitude A. Compute ERSP using newtimef() with 'cycles'=0 and 'scale'='abs'. The integral of powbase over frequency should equal A*A/2, irrespective of changes to analysis parameters such as winsize, sample rate, padratio, etc.

See attached scripts: modified_newtimef.m, test_newtimef.m

FOR THE WAVELET ANALYSIS:

Within newtimef, the wavelet calculations are performed by calling timefreq(). After this, power spectral density is computed. The relevant line of code is 1192:

P  = alltfX.*conj(alltfX); % power for wavelets

I recommend changing this line of code to:

alltfX = sqrt(2) / g.srate * alltfX;

P = real( alltfX .* conj(alltfX) );

% Notes: alltfX returned by timefreq() is a complex array containing
% values PROPORTIONAL to the signal amplitude and its correlation
% with the complex morlet wavelet at each requested time and scale.
%
% 1. Scale alltfX by 1/g.srate. The complex morlet wavelet returned
%    by dipflt3() has not been adjusted for sample rate. Wavelets
%    should (a) sum to zero, and (b) have sum of squares = 1,
%    irrespective of sample rate.
% 2. Scale alltfX by sqrt(2).
%    Not absolutely sure why this is required, but suspect that just
%    as with FFT analysis, the complex morlet wavelet extracts power
%    for both positive and negative frequencies, and the POWER
%    associated with both needs to be summed.
%
% P is now a real array containing the Power Spectral Density of
% the signal. Note that wavelet analysis is NOT identical to frequency 
% analysis, however, when using Morlet wavelets, there is a strong 
% similarity. In this case, P can be considered to have units uV^2/Hz, 
% and the integral (area under the curve) of P equals the total power 
% (variance) of the original signal (providing frequency values are
% sampled adequately).

[arnodelorme]

Test code and modified function from Ross Fulham

Test_Modified_newtimef.zip

[arnodelorme]

Dear Ross, Thanks for your comments. There is actually a long-standing-still-going discussion about this issue on the EEGLAB list. You might be interested in reading it. The solution was assigned to collaborator a while ago. Maybe you might want to contribute with your comments and finally get to the final fix?

See the thread:'[Eeglablist] Units in output of timefreq - wavelet normalization' in this link https://sccn.ucsd.edu/pipermail/eeglablist/2016/thread.html#11562

Thanks for your help on this, Ramon

[arnodelorme]

Ramon,

Would you mind to test the attached function?

Arno

[nucleuscub]

@arnodelorme, the function doesn't brake at least. Now, I see the code proposed here is basically rolling back the method to use the normalization we used before the one currently in use. The old code has been uncommented and the current method deleted. This is a sensitive issue way beyond my expertise ... I suggest to reopen the discussion instead of rolling back to the old method...