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.
```