Finding FWHM of wavelet, and questions on wavelet defs for Figure 18.2

I’m using a mashup of code from ANTS chapters 13 and 18 to analyze my EEG data.

1.

I would like to find the FWHM of the wavelet I’ve created, but cannot find the correct X scaling to use.
The first wavelet is created at 3Hz, 3 cycles, time range is -2 to 2secs, sample rate 1kHz. Length of my EEG single-channel data is 62500 samples, so this makes n_conv_pow2 131072 when length of EEG is added to length of wavelet.
As can be seen in the figure attached, the wavelet is in fact 3Hz and 3cycles, which I also verified with the SignalAnalyzer app in MATLAB. And in the right figure, the peak seems to be at ~98. Somehow, that relates to 3Hz.
My understanding is that the fft will return frequencies from 0 to 131072/2+1, so the 3Hz bin should be at offset 4, rather than 98.

2.

Also, I’m curious about the A scaling factor ( pi * freq * sqrt ( pi ))^- 0.5 that is used, which comes from the example code for Figure 18.2. From what I read in ANTS (pg. 158), this isn’t needed when baselining data.

3.

Also, am I correct that the / freq at the end of the wavelet definition is for 1/f scaling?

Any help much appreciated!

freq_min = 3;
freq_max = 30;
freq_num = 30;
baselinetime = [-800 -200];
channel = 'C3';


channel_ndx = strcmpi(channel, {EEG(1).chanlocs.labels});

frequencies   = logspace(log10(freq_min), log10(freq_max), freq_num);
srate   = EEG.srate;
time    = -2 : 1/srate : 2;                                                 % pg. 168
half_of_wavelet_size = (length(time)-1)/2;

n_wavelet   = length(time);
hz      = linspace(0,srate/2,floor(n_wavelet/2)+1);
n_data      = EEG.pnts;
n_convolution = n_wavelet + n_data - 1;
n_conv_pow2     = pow2(nextpow2(n_convolution));

wavelet_cycles = logspace(log10(3), log10(10), length(frequencies));

fig = figure;
fig.Position=[-2600 200 1600 1200];

baselined_ftt = zeros(length(frequencies), EEG.pnts, EEG.trials);

yplots = 2;
xplots = 2;

for trial = 1 : EEG.trials
    
    % FFT of data
    fft_data = fft(squeeze(EEG.data(channel_ndx, :, trial)), n_conv_pow2);

    % Initialize output time-frequency data
    tf_data = zeros(length(frequencies), EEG.pnts);

    for fi = 1 : length(frequencies)
        
        freq = frequencies(fi);
        cycles = wavelet_cycles(fi);

        % create wavelet and get its FFT  
        wavelet = ( pi * freq * sqrt ( pi ))^- 0.5  * ...                   % scaling factor A see pg. 157
             exp ( 2*1i * pi * freq .* time ) .* ...                        % complex sine wave
             exp (-time .^ 2.0 / ( 2*( cycles  / ( 2 * pi * freq ))^2)) ... % Gaussian
             / freq;                                                        % 1/f scaling
        fft_wavelet = fft(wavelet, n_conv_pow2);                            % pg. 134
        
        % FWHM stuff
        fwave = fft_wavelet;
        fwave = abs(fwave(1:length(hz)))*2;
        fwave = fwave - min(fwave);
        fwave = fwave / max(fwave);
        [~, peakx] = max(fwave);
        [~, left5] = min(abs(fwave(1:peakx)-0.5));
        [~, right5] = min(abs(fwave(peakx:end)-0.5));
        right5 = right5 + peakx - 1;
        
        subplot(yplots, xplots, (fi-1)*1+1);
        plot3(time,real(wavelet), imag(wavelet))
        xlabel('time'), ylabel('real'), zlabel('imaginary')

        subplot(yplots, xplots, (fi-1)*1+2);
        plot(hz, fwave, '.-');
        hold on;
        plot(hz(left5), fwave(left5), 'ro', 'MarkerSize', 10);
        plot(hz(right5), fwave(right5), 'ro', 'MarkerSize', 10);
        plot([hz(left5) hz(left5)], [0 fwave(left5)], 'r');
        plot([hz(right5) hz(right5)], [0 fwave(right5)], 'r');etc.

Hi John. First of all, check out this paper for a deeper discussion of characterizing FWHM.

It looks like you’re doing the right thing, just make sure you are converting the FWHM in indices back into Hz. You can create a vector of frequencies as hz = linspace(0,EEG.srate/2,nConvPnts); and then use that to plot the spectrum. Then the FWHM will be something like diff(hz([right5 left5])).

As for the normalization factors (your questions 2 and 3), please ignore them. They are kindof arbitrary and one of my regrets about that book is not discussing that in more detail. Wavelet normalization in the time domain is really difficult and annoying, whereas normalization in the frequency domain is trivially easy – just divide the spectrum by its maximum value:
fft_wavelet = fft_wavelet / max(abs(fft_wavelet));

1 Like