Introduction
Neurobiologists often want to detect and characterize action potentials bursts. A simple way of looking for the burst period, plus higher frequency details, is to fit the burst structure using Fourier smoothing. A waveform file which represents bursts of action potentials is rectified (positive and negative spikes count equally), then sent to an FFT. The resulting complex spectrum is clipped to low frequencies and inverted back to the time domain. The value of the clipping frequency picks the amount of detail to be included in the result. The time domain version of the band-limited spectrum is then a "best fit" up to the clipping frequency.
The program
The following matlab program implemented the algorithm.
%Spike analysis
%quantify burst duration
%FFT AP burst fit
clear all
clf
dt = 0.0001;
%A file consisting a voltage trace and associated times
load 'C:\My Documents\BigData\bestman\t35-02.mat'
%Choose to analyse from time=41.12 seconds to end of file
v = t35_02(fix(41.1233/dt):end,2);
t = t35_02(fix(41.1233/dt):end,1);
%Center the data around zero
v = v - mean(v);
%Threee plots
np=3
%The data
subplot(np,1,1)
plot(t,v,'color',[0,.75,.75]);
title('Raw spikes')
hold on
%Taking abs value means that + or - spikes add to result
v = abs(v);
%detail parameter: add about 5 for each fourier component that you want
detail = 10;
subplot(np,1,2)
%do the FFT and select just the lowest frequencies according to 'detail'
freqv=fft(v);
freqv(detail:end-detail)=0;
sinefit = ifft(freqv);
plot(t,v,'color',[0,.75,.75]);
hold on
%mult by 3 for clarity on plot
plot(t,real(sinefit)*3,'k','linewidth',2);
title('FFT fit -- first 10 frequencies')
subplot(np,1,3)
%do the FFT and select just the lowest frequencies according to 'detail'
detail = 40;
freqv=fft(v);
freqv(detail:end-detail)=0;
sinefit = ifft(freqv);
plot(t,v,'color',[0,.75,.75]);
hold on
plot(t,real(sinefit)*3,'k','linewidth',2);
title('FFT fit -- first 40 frequencies')
The Result
The top panel is a portion of the input file. The second panel shows an overlay of the rectified data and the first 10 frequency components of the FFT fit. Since the fundamental frequency is about 5 bursts, 10 frequency components corresponds to just two siginficant sine harmonics. From this panel, it is easy to get the average burst period. The third panel shows an overlay of the rectified data and the first 40 frequency components of the FFT fit, or about 8 harmonics of the fundamental burst frequency. The peak shapes are nicely captured.