Differences

This shows you the differences between two versions of the page.

--- analysis:nsb2015:week6 [2015/07/19 21:12]
mvdm created
+++ analysis:nsb2015:week6 [2018/07/07 10:19] (current)
@@ Line 1: / Line 1: @@
 ~~DISCUSSION~~
-===== Filtering: filter design, use, caveats =====
+===== Fourier series, transforms, power spectra =====
 Goals:
-  * Familiarize yourself with basic filtering concepts: frequency and phase response, difference equation, roll-off, ripple
+  * Refresh your memory of basic trigonometry and become comfortable plotting periodic functions in MATLAB
-  * Learn to use MATLAB's filter design and visualization tools
+  * Understand the central idea of the Fourier transform and implement your own decomposition and reconstruction code
-  * Understand the tradeoffs inherent in filtering and use this knowledge to select the appropriate filter for a particular application
+  * Appreciate the fundamental issue of spectral leakage, and how analysis parameters and signal properties affect it
+  * Gain a basic familiarity with different methods for spectral estimation, and use some nonparametric estimators to obtain power spectral densities (PSDs) for some interesting neural signals
-Exercise:
+Resources:
-  * A function ''detectSWR()'' that can reliably identify sharp wave-ripple complexes in a hippocampal LFP using filtering
+  * (bakcground) {{:analysis:course:weeks_ch4_sinusoids.pdf|Chapter 4}} from Weeks, "Digital Signal Processing Using MATLAB & Wavelets"
+  * (background) {{:analysis:course:weeks_ch6_fourierintro.pdf|Chapter 6}} from Weeks, "Digital Signal Processing Using MATLAB & Wavelets
+  * (optional) [[http://www.mathworks.com/help/signal/ug/spectral-analysis.html|MATLAB's guide to spectral analysis]]
-Resources:
+==== Introductory remarks ====
-  * (background) Leis Section 3.10 (difference equation)
+From the very first electrical recordings of the human brain (Berger, 1929) it has been clear that //oscillations// are a salient feature of brain activity. The functional relevance, underlying mechanisms and clinical applications of oscillatory activity continue to be active areas of research; all of these endeavours rely critically on accurate characterizations and quantifications of the data. A fundamental tool in the analysis of oscillations is Fourier analysis, which is introduced in this module.
-  * (background) Leis Chapter 8 (filtering; note that this mentions some concepts that are beyond the scope of this course, so skim over this)
-===== Introductory remarks =====
+==== Generating and plotting basic sinusoids ====
-In the previous module we saw that any signal can be decomposed into a sum of sinusoids, described by a series of magnitude and phase coefficients of a harmonic series. The Fourier transform and the algorithm that performs it can be thought of as mapping the signal from the time domain into the frequency domain. The inverse Fourier transform does the opposite.
+To study oscillations we require //periodic// functions of time, such as ''sin(t)'', which repeat its values in regular intervals.
-This raises the possibility of manipulating the signal in the frequency domain, for instance by removing or amplifying certain frequencies, and then reconstructing the signal. This is an intuitive way to think about //filtering//, defined as an operation or process that removes or attenuates certain features from a signal. [[http://cns-alumni.bu.edu/~slehar/fourier/fourier.html | This page]] has a nice graphical illustration of Fourier filtering.
+Recall that to plot this sort of function in matlab, we first define the time axis (commonly a variable with the name ''tvec''), pass this to the function we wish to plot as an argument, and plot the result:
-Filtering is of central importance in neuroscience, both as an analysis tool (especially when dealing with continuously sampled data) and as a model for operations performed by neural circuits at multiple levels. For instance, the classical receptive fields of V1 neurons can be thought of as filters operating on visual input, and the characteristic time course of postsynaptic potentials imposes limits on how fast signals can be transmitted. Here, we focus on some basic data analysis applications.
+<code matlab>
+%% plot a simple sinusoid
+Fs = 100; % in samples per second (Hz)
+t0 = 0; t1 = 1; % start and end times
+tvec = t0:1./Fs:t1; % construct time axis
-===== Step-by-step =====
+f = 2; % frequency of sine to plot
+y = sin(2*pi*f*tvec); % note sin() expects arguments in radians, not degrees (see sind())
-Filtering is a complex topic in signal processing, with a huge literature, much current research, and many mathematical derivations beyond the scope of this course. However, it is important to be familiar with some of the fundamentals so you can better place descriptions you encounter in the literature, as well as be aware of issues relating to your own analyses. Thus, we will begin with a brief conceptual overview.
+stem(tvec,y);
+</code>
-==== Basic concepts in filtering ====
+Apart from its frequency (2Hz in the above example), a sinusoid is further characterized by its //phase// and //amplitude//, which are readily manipulated as follows:
-=== The difference equation ===
+<code matlab>
-Removing specific frequencies in the Fourier domain is an intuitive way to think about filtering, but this is not how filters are generally implemented. (There are several reasons for this, such as the difficulty in filtering real-time as samples are coming in.) Instead, digital filtering is typically accomplished with a //difference equation// of the form:
+phi = pi/2;
-{{ :analysis:course:diffeq.png |}}
+figure;
+y = sin(2*pi*f*tvec + phi); % a phase shift
+stem(tvec,y);
+hold on;
+plot(tvec,cos(2*pi*f*tvec),'LineWidth',2) % notice, cosine is simply phase shifted sine
+legend('Sin (Phase Shifted)', 'Cos');
+hold off;
-This equation describes how to compute sample //n// of the filtered signal //y// from the original signal //x//. The value of //y(n)// in general can be a function of an arbitrary number of past samples (-1, -2, etc.) of the original signal //x//, as well as of //y// itself.
+a = 2;
-The coefficients ''b'' and ''a'' fully define the filter. Notice that a filter with only ''b'' coefficients depends only on the original signal //x//; that is, there are no //feedback// terms taken from the filtered signal //y//. The ''a'' coefficients describe the feedback components. This distinction is the basis for the commonly used terms "FIR" (Finite Impulse Response; only ''b'' components) and "IIR" (Infinite Impulse Response, ''a'' components) to describe a filter. The behavior of IIR filters can be much more complex because of this feedback.
+figure;
+y = a.*sin(2*pi*f*tvec + phi); % amplitude change
+stem(tvec,y);
-Another useful term to be aware of is the //order// of a filter: this simply refers to the maximum number of samples back it looks. This is equivalent to the maximum number of coefficients ''a'' and ''b'' the filter has; ''a(0)'' and ''b(0)'' always need to be defined (''a(0)'' is implicit on the left side in front of ''y(n)'', often omitted because it is generally 1) as well as at least one ''a'' or ''b'' for each sample back up to the maximum.
+</code>
-(Technical note: the difference equation can be related to the frequency domain through the Z-transform; it is not necessary to understand this, but if you are curious, the Leis book has an explanation.)
+As an aside, sinusoids are the basis for wireless transmission, such as that for FM radio, where a message is encoded in the frequency modulation of a "carrier" signal:
-=== An example using filter() ===
+<code matlab>
+%% frequency modulation (FM) signal
+f2 = 10;
+m = 2;
-An example of a simple filtering operation is to compute a running average. Looking at the difference equation above, we see we can accomplish this using ''b'' coefficients only. How many coefficients we use will determine the size of the window (in samples) that we average over; the magnitude of the coefficients should be set so that we in fact get the mean, and not for instance the sum.
+subplot(311)
+s1 = sin(2*pi*f*tvec);
+plot(tvec,s1); title('message');
-So, to compute the running mean over four samples, we want:
+subplot(312);
+s2 = sin(2*pi*f2*tvec);
+plot(tvec,s2); title('carrier');
-{{ http://www.mathworks.com/help/releases/R2013b/matlab/data_analysis/eqn1133878950.png?300 }}
+subplot(313);
+s3 = sin(2*pi*f2*tvec + m.*sin(2*pi*f*tvec - pi/2));
+plot(tvec,s3); title('FM signal');
+</code>
-Using ''filter()'', we would simply do
+Notice how the sinusoid in the bottom plot smoothly modulates its frequency. Many oscillations in the brain, such as hippocampal theta and ventral striatal gamma, show such modulations as well.
-<code matlab>
+==== Sums of sinusoids and harmonic series ====
-a = 1; % a_0 is the (hidden) coefficient on the left side, in front of y(n)
-b = [1/4 1/4 1/4 1/4]; % four b's of 1/4 each so we get the mean
-y = filter(b,a,x); % x is the original signal, y the filtered version
+What happens if we sum some sinusoids together? A compact representation for this is a //harmonic series//:
-</code>
-Of course, this won't work because we don't have our input signal ''x'' defined yet. As a quick illustration (following [[http://www.mathworks.com/help/matlab/data_analysis/filtering-data.html | this]] MATLAB page), we can do:
+$$S(t)=\sum_{n=1}^{N}{a_{n} cos(2 \pi n f_{0} t+\phi_{n})}$$
+Note that each sinusoid of this series has a frequency that is an integer multiple $n$ of some base frequency $f_0$. This is a compact representation because we only specify the amplitude $a_n$ and phase $\phi_n$ for each of the series, with the frequencies fixed.
+Figure 4.7 in the Weeks chapter is a sum of four sinusoids, with a base frequency of 2:
 <code matlab>
-load count.dat;
+%%Harmonic series example
-x = count(:,1);
+mag = [0.1 0 1.3 0.5]; % magnitudes for each term
+pha = [-pi/6 0 pi 2*pi/3]; % phases for each term
+f = 2; % base frequency
-t = 1:length(x);
+signal_out = zeros(size(tvec));
-plot(t,x,'-.',t,y,'-'), grid on
+for ii = 1:numel(mag) % note, the book chapter uses i, not best practice!
-legend('Original','Filtered',2)
+    this_signal = mag(ii)*cos(2*pi*f*ii*tvec + pha(ii));
+    plot(tvec,this_signal,'r:'); hold on;
+    signal_out = signal_out + this_signal; % build the sum
+end
+figure;
+plot(tvec,signal_out,'LineWidth',2);
 </code>
-Combining the above pieces should give:
+☛ Why is it not a good idea to use the variable name ''i''? (Hint: for the same reason you should not define a variable named ''pi''...)
-{{ :analysis:course:week5_fig1.png?600 |}}
+It looks like we can create some interesting signals by summing sinusoids.
-The filtered signal looks roughly as expected -- a nicely smoothed version of the original -- but a few things are worth noting. Look at the first sample of the filtered signal: it equals ''x(1)/4'', which means that any values of ''x'' which we did not have were assumed to be zero. This is an example of an "edge effect"; in general we don't really think the signal we didn't sample is truly zero, but this is ''filter()'''s implicit, default assumption.
+In fact, the central insight underlying Fourier analysis is that we can use a sum (series) of sinusoids to approximate //any// signal to arbitrary precision. An important corollary to this is that any signal can be //decomposed// into a series of sinusoids. Let's demonstrate this.
-Another property of this filtered signal is that it is //phase-shifted// to the right; this of course arises because our ''y(n)'' is based only on past samples, not on the future. This is a key issue for neuroscience data analysis and we will return to it below.
+==== Decomposing and reconstructing a signal ====
-=== Common filters and their applications ===
+We will use the MATLAB function ''fft()'' (for **F**ast **F**ourier **T**ransform; a technical point discussed in the Weeks chapter is that this in fact computes something called the Discrete Fourier Transform or DFT) to obtain the amplitudes and phases ($a_n$ and $\phi_n$ in the equation above) of a randomly generated signal, and then piece the signal back together by plugging them in to a harmonic series:
-Computing a running mean has its uses, but for neural data we are typically interested in other applications of filtering. The canonical filter types, illustrated below, are the //lowpass// filter (pass low frequencies, suppress high frequencies), //highpass// filter (the reverse), //bandpass// filter (only pass frequencies within a certain range) and //notch// filter (suppress frequencies within a narrow range):
+<code matlab>
+x = round(rand(1,8)*10); % generate a length 8 vector of integers between 0 and 10
+xlen = length(x);
-{{ http://images.books24x7.com/bookimages/id_15261/p30001dd3g4320001vpp_thm.jpg?nolink&600 |}}
+% get magnitudes and phases of Fourier series
+X = fft(x);
+Xmag = abs(X); % magnitudes, a_n
+Xphase = angle(X); % phases, phi_n
-Note that these illustrations are in the frequency domain: some frequencies are unchanged (0 dB; recall that the decibel is a common unit of signal power), while others are attenuated. The -3dB point is a common reference and corresponds to a 50% reduction in signal power (recall that dB is a log scale).
+n = 0:xlen-1;
+t = 0:0.05:xlen-1; % a finer timescale to show the smooth signal later
-How do we know what values of ''b'' and ''a'' will accomplish these filtering operations? Can the //optimal// ''b'' and ''a'' for a given application be found? These are questions about //filter design//, a field for which MATLAB (in particular the Signal Processing Toolbox) provides many useful tools.
+for iH = xlen-1:-1:0 % reconstruct each harmonic
+    s(iH+1,:) = Xmag(iH+1)*cos(2*pi*n*iH/xlen + Xphase(iH+1))/xlen;
+    sm(iH+1,:) = Xmag(iH+1)*cos(2*pi*t*iH/xlen + Xphase(iH+1))/xlen;
+    % detail: xlen appears here because the fundamental frequency used by fft() depends on this
+end
-As may be expected from your experience with windowing and spectral leakage in the previous module, there is no ideal filter that completely rejects all unwanted frequencies while leaving all desired frequencies intact. Rather, different filters make different tradeoffs in terms of their frequency and phase response.
+ssum = sum(s);
+smsum = sum(sm);
-For instance, a //Butterworth// filter has a maximally flat frequency response in the passband. However, its //rolloff// (the transition from passband to stopband) is not as steep as that of //Chebyshev// filters. In turn, Chebychev filters experience //ripple// (some distortion) in either the passband or the stopband and thus are not as flat as Butterworth filters. Many other filters, with different properties, exist. Thus, it becomes important to select the right filter for your application. To help with this process, MATLAB has a useful tool that provides the key properties of a given filter at a glance.
+figure;
+plot(n, x, 'go', t, smsum, 'b', n, ssum, 'r*');
+legend({'original','sum - all','sum - points only'});
+</code>
-=== Frequency response of some common filter types ===
+This gives something like:
-Let's start with designing a basic Butterworth bandpass filter. The ''help'' for ''butter()'' says:
+{{ :analysis:course:week4_fig1.png?nocache&600 }}
-<code matlab>
+Notice that the reconstruction is essentially perfect. We retrieved the original 8 points by plugging in the coefficients returned by ''fft()'' into a series of sinusoids.
->> help butter
- butter Butterworth digital and analog filter design.
-   [B,A] = butter(N,Wn) designs an Nth order lowpass digital
-   Butterworth filter and returns the filter coefficients in length
-   N+1 vectors B (numerator) and A (denominator). The coefficients
-   are listed in descending powers of z. The cutoff frequency
-   Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to
-   half the sample rate.
-   If Wn is a two-element vector, Wn = [W1 W2], butter returns an
-   order 2N bandpass filter with passband  W1 < W < W2.
-</code>
-So, ''butter()'' computes for us the coefficients ''a'' and ''b'' of the difference equation. The "descending powers of z" referred to correspond to samples further in the past: z^0 is the current sample, z^-1 the previous one, et cetera. This is Z-transform stuff you don't need to understand; for now it's sufficient to know that ''butter()'' returns the ''a'' and ''b'' coefficients in the order we expect. Likewise, the numerator and denominator refer to where ''a'' and ''b'' would end up if we wrote out the "transfer function" for the filter; no need to worry about this either.
-Note that we need to specify the cutoff frequencies as numbers between 0 and 1, where 1 correponds to the Nyquist frequency of our data. So we cannot directly say we want the cutoff to be e.g. 250Hz, we have to normalize by ''Fs/2''.
+☛ What is the purpose of running the ''for'' loop backwards in the above code? (Hint: recall that MATLAB slows down substantially if it needs to allocate memory for expanding variables.)
-☛ Generate a 10-second long white noise signal, sampled at 500Hz. Filter it using a bandpass Butterworth filter between 50 and 100 Hz of order 4. Plot the Welch spectrum, in dB, of the original and filtered signal, using a 512-sample Hanning window. Evaluate the FFT over 2^14 points.
+The set of magnitudes and phases that describe the harmonic series that reconstruct the signal are known as the //magnitude spectrum// and //phase spectrum// respectively. The magnitude spectrum can be interpreted as signal //power// at different frequencies.
-Your code should look something like:
+==== Interpreting the output of MATLAB's fft() function ====
+You might have noticed the magnitude and phase spectra above have a peculiar form. Let's explore this a little further using a signal with known frequency and phase content:
 <code matlab>
+Fs = 20; % in samples per second (Hz)
+t0 = 0; t1 = 1; % start and end times
+tvec = t0:1/Fs:t1-(1/Fs); % construct time axis; generate exactly 20 samples
-% set up time axis
+f = 2; % signal frequency
-Fs = ...
+y = sin(2*pi*f*tvec); % construct signal, a 2Hz sine wave sampled at 20Hz for 1s
-tvec = ...
-% generate white noise
+yfft = fft(y,length(y));
-x = rand(...)
+yfft_mag = abs(yfft); yfft_ph = angle(yfft);
+stem(yfft_mag)
+</code>
-% get PSD
+The result:
-[Porig,Forig] = pwelch(x, ...)
-% design filter
+{{ :analysis:course:week4_fig2.png?600 |}}
-W1 = ...
-W2 = ...
-[b,a] = butter(4,[W1 W2]);
-y = filter(...)
-% get PSD
+Some issues are apparent:
-[Pfilt,Ffilt] = pwelch(y, ...)
-% plot the resulting PSDs
+  * ''fft()'' did not return a frequency axis, so the output is not straightforward to interpret.
-subplot(121)
+  * Whatever the coefficients returned, there is not a single peak as we would expect from the pure input signal: in fact there are two peaks.
-plot(... 10*log10(..));
+To understand this, recall that the largest frequency that can be detected in a signal sampled at ''Fs'' is ''Fs/2'', the //Nyquist frequency//. Thus, we would expect the frequency axis from ''fft()'' to go from 0 to ''Fs/2'' at most. In addition, it turns out the Fourier transform is defined for //real// as well as for //complex// signals, and it returns spectra for both these components. Since we are only interested in real signals, we get a symmetric spectrum back (note that if we did not use ''abs()'' on the output of ''fft()'' we would get some imaginary components as well).
+To construct a frequency axis that takes both these ideas into account, we can do:
+<code matlab>
+Npoints = length(y);
+F = [-Npoints/2:Npoints/2-1]./Npoints; % construct frequency axis
+yfft_mag = fftshift(yfft_mag); % align output, see note below
+stem(F,yfft_mag);
+xlabel('Frequency (Fs^{-1})');
 </code>
-When done correctly, the resulting PSDs should be similar to these:
+''fftshift()'' cuts the second (complex) half of the spectrum and pastes it backwards at the beginning, so that our frequency axis is now correct; it is in units of ''1 / Fs'' so 0.1 corresponds to the 2Hz we put in.
+For analysis purposes we don't care about the complex part of the spectrum, and this is usually ignored to yield what is referred to as the "single-sided spectrum". Because we would never actually use ''fft()'' directly to estimate spectral content (see the section on spectral estimation below) we won't do more work to get rid of it now.
+Notice also the superscript in the ''xlabel'': MATLAB can interpret basic [[http://web.ift.uib.no/Teori/KURS/WRK/TeX/symALL.html|LaTeX math symbols]] for figure text.
-{{ :analysis:course:week5_fig2.png?600 |}}
+A final point about the output of ''fft()'' is that it contains the magnitude spectrum at specific frequencies. For instance, we have the 0.1*Fs point, but not 0.125*Fs! This will become important later.
-As you can see, our filter is doing something, but it's also clear that it's not perfect. Frequencies outside the passband still get passed to some degree, and if you look carefully (''grid on'' can help), you can see that frequencies in the passband but close to the rolloff frequencies are slightly attenuated already.
+==== Zero-padding the FFT ====
-Of course, in general we are not interested in white noise, but it is a useful testbed to gauge the properties of a filter. Because we know white noise has a flat frequency spectrum, we can see at a glance which frequencies are attenuated after filtering.
+The above example was constructed nicely so that our signal contained exactly two full periods. This will not be true for real world signals. What happens if we don't have this perfect input?
-One way to improve is to ask MATLAB to suggest an appropriate filter order for us:
+☛ Change the ''tvec'' variable above to contain one more sample, like so:
 <code matlab>
-Wp = [ 50 100] * 2 / Fs; % passband - between 50 and 100 Hz
+tvec = t0:1/Fs:t1;
-Ws = [ 45 105] * 2 / Fs; % stopband
-[N,Wn] = buttord( Wp, Ws, 3, 20); % determine filter parameters
-[b2,a2] = butter(N,Wn); % builds filter
 </code>
-The ''buttord()'' function takes the filter specifications and returns a suggested filter order (''N'') and new frequency cutoffs ''Wn'' to feed to ''butter()''. The way we specify this is to say that we require a minimum level of attenuation in the stopband (in this case, 20dB) and we are willing to tolerate a certain amount of distortion ("ripple") in the passband (in this case, 3dB).
+This means that our signal is now no longer an integer number of periods. The resulting two-sided spectrum is:
-As it turns out, for this filter, ''buttord()'' suggests order 15! This is quite a difference from the 4th order filter we implemented above. Let's see how our new filter compares. Happily, we don't need to keep filtering white noise in order to see the frequency response of a filter, because MATLAB provides a nice tool:
+{{ :analysis:course:week4_fig3.png?600 |}}
+Notice that:
+  - The peaks now appear at a frequency not exactly equal to the true frequency
+  - Other non-zero components have appeared
+To explain (1), recall that ''fft()'' evaluates the magnitudes of specific frequencies only. Inspection of the spectrum above indicates that we don't have a frequency bin that is exactly 0.1 (the true frequency of our signal). Let's fix that using the second, optional, argument of ''fft()''.
 <code matlab>
-fvtool(b,a,b2,a2)
+tvec = t0:1/Fs:t1;
-</code>
+nPoints = [length(tvec) 64 256 1024];
-You should get:
+for iP = 1:length(nPoints) % repeat fft with different numbers of points
-{{ :analysis:course:week5_fig3.png?600 |}}
+    nP = nPoints(iP);
+    subplot(2,2,iP);
+    y = sin(2*pi*f*tvec);
+    yfft = fft(y,nP);
+    yfft_mag = abs(yfft); yfft_ph = angle(yfft);
+    F = [-nP/2:nP/2-1]./nP;
+    yfft_mag = fftshift(yfft_mag);
+    plot(F,yfft_mag,'kx',F,yfft_mag,'k');
+    title(sprintf('%d point FFT',nP));
+    xlabel('Frequency (Fs^{-1})');
-Notice how our new filter (in green) is much more effective than the previous one. It has sharper roll-off and better attenuation in the stopband. The units on the frequency axis are fractions of ''Fs/2'', so 0.2 corresponds to 50Hz as expected.
+end
+</code>
-☛ What happens if you get greedy and try to have a stopband of [48 102] as an input to ''buttord()''?
+This gives:
-Let's try a different filter, a Chebyshev Type I. With this one, we can be greedy:
+{{ :analysis:course:week4_fig4.png?600 |}}
-<code matlab>
+As we increase the number of points to evaluate the FFT over, we get increased frequency resolution, and indeed the peaks of the spectrum converge to 0.1 as expected. Under the hood, this is in fact accomplished by //padding the input signal with zeros// before the DFT is computed; doing this does not change the spectral content (the zeros are not periodic), but allows a much longer harmonic series, with a smaller fundamental frequency, to be evaluated.
-Wp = [ 50 100] * 2 / Fs;
-Ws = [ 48 102] * 2 / Fs;
-[N,Wn] = cheb1ord( Wp, Ws, 3, 20);
-[b_c1,a_c1] = cheby1(N,0.5,Wn);
-fvtool(b2,a2,b_c1,a_c1)
-</code>
-Note that we use the same workflow of having MATLAB suggest a filter order and passband based on our specifications. The ''cheby1()'' function needs one additional input argument compared to ''butter()''; this relates to the "ripple" that is visible in the frequency response:
+A typical value to use for the number of points to evaluate the FFT is the next power of 2 (after however many samples your signal contains). This is because part of what makes the FFT fast is that it can easily divide the signal in half. But non-power of 2 values also work.
-{{ :analysis:course:week5_fig4.png?600 |}}
+☛ What happens if you try to evaluate the FFT using a number of points //smaller// than that in your signal?
-As you can see, our Chebyshev filter (in green) has a sharper rolloff, but at a cost: the frequency response in our passband is now no longer flat, and has a "ripple" instead. There is also a Chebyshev Type II filter, which has a flat passband response but a ripple in the stopband, but its rolloff tends to be less sharp so is less commonly used.
+As we increase the number of points to evaluate the FFT, we can obtain coefficients for frequencies of arbitrary precision. For this reason, the power spectrum is often referred to as //power spectral density// or PSD.
-=== Phase responses and filtfilt() ===
+FIXME Demonstrate that you get the same effect by zero-padding (rather than making the signal longer)
-Let's apply our filter to a more realistic signal:
+==== Spectral leakage ====
+It is clear from the above figure that even if we increase the frequency resolution of our FFT with zero-padding, we still have an imperfect estimate of the true frequency content of our signal. In particular, the estimate around the true frequency has a nonzero width -- this part of the magnitude spectrum is referred to as the //main lobe// -- and we have nonzero spectral content for other frequencies as well (the //side lobes//).
+To explore the source of this **spectral leakage**, let's make our signal longer while evaluating the FFT over the same number of points:
 <code matlab>
-Fs = 500; dt = 1./Fs;
+tvec = t0:1/Fs:t1-(1/Fs);
-t = [0 10];
+nRepeats = [1 2 4 8];
-tvec = t(1):dt:t(2)-dt;
-s1 = sin(2*pi*80*tvec+pi/6);
+nP =  1024;
-s2 = sin(2*pi*40*tvec);
-s = s1 + s2;
-sf = filter(b_c1,a_c1,s);
+for iP = 1:length(nRepeats)
+    subplot(2,2,iP);
+    y = sin(2*pi*f*tvec);
+    y = repmat(y,[1 nRepeats(iP)]); % repeat the signal a number of times
+    yfft = fft(y,nP);
+    yfft_mag = abs(yfft); yfft_ph = angle(yfft);
+    F = [-nP/2:nP/2-1]./nP;
+    yfft_mag = fftshift(yfft_mag);
+    plot(F,yfft_mag,'kx',F,yfft_mag,'k');
+    title(sprintf('%d repeats',nRepeats(iP)));
+    xlabel('Frequency (Fs^{-1})');
-plot(tvec,s,'k',tvec,sf,'r--'); hold on;
+end
-legend({'original','filtered'});
-xlim([0 0.2]);
 </code>
 The result:
-{{ :analysis:course:week5_fig5.png?600 |}}
+{{ :analysis:course:week4_fig5.png?600 |}}
-The filter was effective in removing the lower-frequency (40Hz) component, with only the 80Hz oscillation remaining. However, the //phase// of the signal has clearly changed also, by what appears like 180 degrees -- the faster-oscillation peaks in the original trace now are closely aligned with the troughs in the filtered trace. As was the case in our first filtering example (the moving average filter above), the filtered signal appears delayed relative to the original.
+Notice that the magnitude spectrum converges to the true frequency as we increase the length of the signal. However, for any signal of finite length, there will always be some spectral leakage. This occurs because we are effectively cutting off the signal abruptly; other than changing the length of the signal (which is not always an option) we can attempt to minimize spectral leakage through using less abrupt cutoffs, discussed in the next section.
-Clearly, such phase shifts can be devastating for the analysis of neural data. If features of a LFP are delayed because of filtering, this may obscure relationships between the LFP and behavioral or neural events. In addition, any analysis that relies on knowing the phase of a LFP, such as theta phase precession, or cross-frequency coupling, will be affected as well. In general, the //phase response// is an important characteristic of any filter, and indeed ''fvtool'' can display it.
+==== Windowing ====
-☛ Run ''fvtool'' again on the Butterworth and Chebyshev filters above, and now select the Phase Response button in the top left of the window.
+Computing a FFT over a window of finite size is as if we are taking a finite-size window (for instance '' w(n) = 1 if 0 <= n <= N; 0 otherwise''; note that this defines a //rectangular// window) and multiplying this with a hypothetical infinite signal.
-Notice that the phase response is not constant but in fact depends on the input frequency. This makes it very difficult to correct for such phase shifts. For neural data, where even small phase shifts can be problematic, we therefore take an alternative approach: we filter the signal forwards and backwards, such that the net phase response is zero: no phase shift! This is accomplished using the ''filtfilt()'' function:
+It turns out that the spectrum of the windowed signal equals the //convolution// of the signal's spectrum and the window's spectrum; this occurs because of the [[http://en.wikipedia.org/wiki/Convolution_theorem|convolution theorem]] that states that multiplication in the time domain equals convolution in the frequency domain.
+Thus, it becomes important to understand the spectral properties of the windows used for Fourier analysis. Let's plot a few:
 <code matlab>
-sf = filtfilt(b_c1,a_c1,s);
+nP = 25;
+nPFFT = 1024;
-plot(tvec,s,'k',tvec,sf,'r--'); hold on;
+windows = {'rectwin','triang','hamming','hanning','blackman'};
-legend({'original','filtered'});
+cols = 'rgbcmyk';
-xlim([0 0.2]);
-</code>
-☛ Verify that there is no longer any detectable phase shift.
+for iW = 1:length(windows)
+    eval(cat(2,'wn = ',windows{iW},'(nP);')); % make sure you understand this
+    wn = wn./sum(wn);
+    subplot(211); % plot the window
+    plot(wn,cols(iW),'LineWidth',2); hold on;
+    subplot(212);
+    yfft = fft(wn,nPFFT);
+    yfft_mag = abs(yfft); yfft_ph = angle(yfft);
+    F = [-nPFFT/2:nPFFT/2-1]./nPFFT;
+    yfft_mag = fftshift(yfft_mag);
+    h(iW) = plot(F,yfft_mag,cols(iW),'LineWidth',2); hold on;
+end
-So far so good, but what are the consequences for the frequency response of doing this?
+xlabel('Frequency (Fs^{-1})');
+legend(h,windows);
+</code>
-<code matlab>
+☛ What does the ''eval()'' statement in the above code do?
-%% compare freq responses
-Fs = 500; dt = 1./Fs;
-t = [0 10];
-tvec = t(1):dt:t(2)-dt;
-x = rand(size(tvec)); % white noise input
+The result:
-[P,F] = pwelch(x,hanning(512),256,2^14,Fs);
-y1 = filter(b_c1,a_c1,x);
+{{ :analysis:course:week4_fig6.png?600 }}
-[P1,F1] = pwelch(y1,hanning(512),256,2^14,Fs);
-y2 = filtfilt(b_c1,a_c1,x);
+Inspection of the magnitude spectra of the different windows shows that although the rectangular window has the narrowest mainlobe (a desirable characteristic; narrower means a better estimate of the true underlying frequency) it also has the largest sidelobes (an undesirable characteristic; these sidelobe frequencies are not actually present in the signal!). As you can see, different windows represent different tradeoffs between these properties.
-[P2,F2] = pwelch(y2,hanning(512),256,2^14,Fs);
-plot(F,10*log10(P),F,10*log10(P1),F,10*log10(P2));
+☛ Change the y axis to log scale to better visualize the differences (recall that you can use ''get(gca)'' to access the properties of the current axes).
-legend({'original','filter','filtfilt'});
-</code>
-This gives:
+☛ Modify the code in the Spectral Leakage section, above, to use a Hamming window instead of a rectangular window. Verify that the sidelobes are now smaller, at the expense of a slightly wider mainlobe.
-{{ :analysis:course:week5_fig6.png?600 |}}
+FIXME Clarify if the integral of the windows should be 1 (to preserve power estimates) and where this gets done in the code.
+==== Robust spectral estimation methods =====
-As is often the case, the output from ''filtfilt()'' actually has a steeper rolloff than that from ''filter()''. This is because we are effectively filtering twice, an effect that can be approximated by increasing order of the filter (if you were to filter it only once). ''filtfilt()'' tends to be more robust, but it is always a good idea to check your filter on white noise if you have not used it before.
+The above method of spectral estimation, obtaining the Fourier coefficients by applying the DFT to the entire data set, is also referred to as constructing a //periodogram//. MATLAB has a function for this:
-==== Some typical neuroscience applications ====
+<code matlab>
+[Pxx,F] = periodogram(y,[],nP,Fs);
+plot(F,Pxx); xlabel('Frequency (Hz)');
+</code>
-=== Removing 60Hz line noise in the data ===
+Note that this plots the //one-sided// spectrum, and that the units on the frequency axis are now in Hz. In addition, rather than plotting the magnitudes of the sinusoid components as before, we now have //power//, which is simply the magnitude squared.
-Let's try to design a //notch// filter to remove 60Hz line noise, using the familiar method:
+It is easy to change the window:
 <code matlab>
-[b,a] = butter(10, [59 61] * 2 / Fs, 'stop');
+hold on;
-fvtool(b,a);
+[Pxx,F] = periodogram(y,hanning(length(y)),nP,Fs);
+plot(F,Pxx,'r');
 </code>
-As you can see, the frequency response doesn't look good.
+Note again the trading off of mainlobe width against sidelobe magnitude.
-☛ Try some different filter orders and see if you can get the desired notch shape (i.e. attenuation at 60Hz, no attenuation everywhere else).
+It turns out that the periodogram, although usable, is not a very good spectral estimator in some ways; it is //biased//: its variance does not tend to zero as data length goes to infinity ([[http://www.mathworks.com/help/signal/ug/spectral-analysis.html|details]]). This makes later statistical comparisons difficult, and the spectrum look noisy, as we'll demonstrate in a later section.
-This issue is similar to what we encountered when trying to get a sharper rolloff for our bandpass filter. In that case we fixed things by going to a Chebyshev filter. Another method is the following:
+A useful approach to address these issues is to cut the signal up into smaller segments, estimate the spectrum for each segment, and combine the estimates, a process known as [[http://en.wikipedia.org/wiki/Bartlett's_method|Bartlett's method]]. An example of a spectral estimator that uses this approach is //Welch's method//, which uses segments (or windows) that can overlap.
+Type ''help pwelch'' and read the overall description of what this function does. Based on the discussion above, this should make sense.
+Here is an illustration of how ''pwelch()'' performs:
 <code matlab>
-[z,p,k] = butter(10, [59 61] * 2 / Fs, 'stop'); % note, we ask for 3 outputs instead of 2
+Fs = 20; % in samples per second (Hz)
-[sos,g] = zp2sos(z,p,k); % convert to SOS format
+t0 = 0; t1 = 1;
-h = dfilt.df2sos(sos,g); % create filter object
+f = 2;
-fvtool(h);
+nRepeats = 4;
-</code>
-Now we have a good looking notch filter. This so-called "second-order section" format is more numerically precise than the standard difference equation ''[b,a]'' format. When dealing with higher order filters it can make a difference!
+tvec = t0:1/Fs:t1-(1/Fs);
-☛ Test this nice notch filter on white noise using ''filtfilt()''.
+nP =  1024;
+y = sin(2*pi*f*tvec);
+y = repmat(y,[1 nRepeats]);
-=== Detecting movement artifacts ===
+[Pxx,F] = periodogram(y,rectwin(length(y)),nP,Fs);
+plot(F,Pxx);
-Movement artifacts arising from EMG (electrical activity generated by muscles) are common when recording neural signals from behaving subjects. Ideally, these are removed by correct referencing, but this is not always possible. Chewing artifacts, when a rat consumes food pellets, can be particularly pernicious; as can eyeblinks when recording scalp EEG.
+hold on;
+wSize = 40;
+[Pxx,F] = pwelch(y,rectwin(wSize),wSize/2,nP,Fs);
+plot(F,Pxx,'r'); xlabel('Frequency (Hz)');
+</code>
-As an example:
+Note that the Welch spectrum sacrifices some frequency resolution, because of the smaller window that results from cutting up the data, but our estimate is now more robust.
+==== Pitfalls for real-world signals ====
+We are almost ready to apply these methods to some real data. However, we first need to be aware of an issue that often comes up with real data. ''fft()'' and tools that rely on it, such as the spectral estimators in the previous section, assume that the data is sampled at evenly spaced intervals. However, this is only approximately true for Neuralynx data: chunks of 512 samples are sampled at a regular interval (''Fs = 2000'') but the interval between chunks is slightly larger than within, such that the overall Fs is slightly smaller than 2000. Because this difference is tiny, it will not affect our spectral estimate much. However if the difference becomes large our spectrum will be wrong. For instance:
 <code matlab>
-%% cd to R016-2012-10-08 folder first
+Fs = 20; % in samples per second (Hz)
-cfg = [];
+t0 = 0; t1 = 1; f = 2;
-cfg.fc = {'R016-2012-10-08-CSC02b.ncs'};
+nP =  1024;
-csc = LoadCSC(cfg);
+gaps = [5 10 15]; % idx of samples to be removed
-cscR = restrict(csc,1270,1272);
+tvec = t0:1/Fs:t1;%-(1/Fs);
-plot(cscR.tvec,cscR.data)
+y = sin(2*pi*f*tvec);
-</code>
-This is a piece of LFP recorded as this rat was eating. Note the characteristic rhythmic pattern of high-frequency oscillation events that occur approximately 4 times each second:
+subplot(211)
+plot(tvec,y,'k*'); hold on;
-{{ :analysis:course:week5_fig7.png?600 |}}
+yfft = fft(y,nP);
+yfft_mag = abs(yfft); yfft_ph = angle(yfft);
-We could try to //remove// these events from the signal and then pretend they were never there in further analysis, but a more conservative approach is simply to detect them and store the corresponding times so that we can exclude them from subsequent analysis.
+F = [-nP/2:nP/2-1]./nP;
+yfft_mag = fftshift(yfft_mag);
-Let's take a guess at a frequency band that may be able to detect these events:
+subplot(212);
+plot(F,yfft_mag,'k-x'); hold on;
-<code matlab>
+xlabel('Frequency (Fs^{-1})');
-Fs = cscR.cfg.hdr{1}.SamplingFrequency;
-Wp = [ 180 220] * 2 / Fs;
-Ws = [ 178 222] * 2 / Fs;
-[N,Wn] = cheb1ord( Wp, Ws, 3, 20); % determine filter parameters
-[b_c1,a_c1] = cheby1(N,0.5,Wn); % builds filter
-%fvtool(b_c1,a_c1); % remember to check your filter!
+% signal with gaps
+y = sin(2*pi*f*tvec);
+y2 = y;
+y2(gaps) = []; tvec(gaps) = []; % remove
-y = filtfilt(b_c1,a_c1,cscR.data);
+subplot(211);
-plot(cscR.tvec,cscR.data,'b',cscR.tvec,y,'r');
+plot(tvec,y2,'bo'); hold on;
-</code>
-It looks like this filter is picking up something from the chewing events. But we are not interested in a chewing-band signal per se; we want to use it to detect the presence of chewing events. Thus, we can convert the oscillating signal into an unsigned quantity, signal //power//:
+yfft = fft(y2,nP);
+yfft_mag = abs(yfft); yfft_ph = angle(yfft);
-<code matlab>
+F = [-nP/2:nP/2-1]./nP;
-chew_power = y.^2;
+yfft_mag = fftshift(yfft_mag);
+subplot(212);
+plot(F,yfft_mag,'b-x');
+legend('y','y2 (with gaps)')
 </code>
-Plot this //instantaneous signal power// and notice this is a pretty variable quantity. This is where our moving average filter comes in handy, but we can also use ''medfilt1()'', a median filter which is a bit more robust to outliers:
+Notice that the estimate for the data with some samples removed (in blue) is substantially lower than the true frequency.
-<code matlab>
+In general, if you have unevenly spaced samples, you can consider:
-chew_power = y.^2;
+  * Interpolating the data on an evenly spaced timebase using a function such as ''interp1()''. This is only safe if the frequencies you will be analyzing are of a much lower order than the variation in the sample spacing -- the interpolation can introduce high-frequency components into the signal.
-chew_power_filtered = medfilt1(chew_power,101); % filter window is specified in samples, so this is ~50ms
+  * Using spectral estimators designed to work with uneven spacing, such as the [[http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m|Lomb-Scargle periodogram]].
-[h1 h2] = plotyy(cscR.tvec,cscR.data,cscR.tvec,chew_power_filtered);
-</code>
-The resulting //envelope// is well-behaved and matches the chewing events nicely:
+An extension of the uneven sampling problem is gaps in the data. This can occur if there is a technical problem with the recording system, if there are artifacts that have been removed, or simply if as part of your setup you naturally have gaps (e.g. intertrial intervals). Of course you may decide to cut up your data and analyze each segment (with even sampling) separately, much like Bartlett's method itself. However often it is convenient not to have to do this.
-{{ :analysis:course:week5_fig8.png?600 |}}
+A reasonable solution to deal with gaps is to fill them with zeros, such that a constant ''Fs'' is maintained. The zeros will introduce edge effects due to windowing. In addition, if there is temporal structure in where the gaps appear this will also show up in the spectrum: for instance, imagine there was a repeating pattern of 100ms gap, followed by 100ms of data, then a strong 10Hz component would appear. Thus, if you are filling gaps in this way, it is prudent to check the //autocorrelation// of the times where gaps occur, to form an idea of where artifactual components might show up; we will treat this later in the course.
+==== Application to real data ====
-FIXME update image file above to reflect use of plotyy()
+Let's begin by loading a ventral striatal LFP signal:
-The filtered chewband power can now be used for a simple thresholding operation to decide if a chewing event is present or not. By changing the filter properties, the same approach can be used to obtain a time series of any frequency band of interest, to examine for instance the relationship between running speed and theta power, or between gamma events and reward anticipation.
+<code matlab>
+% cd to your R016-2012-10-08 folder
+cfg = [];
+cfg.fc = {'R016-2012-10-08-CSC04d.ncs'};
+csc = LoadCSC(cfg);
+</code>
-This approach illustrates that often we don't really want to //replace// the original signal with a filtered version. Rather, the same original signal is used to generate multiple different filtered signals, used together or separately to explore distinct analysis questions and to make include/exclude decisions.
+This dataset contains a baseline (rest) recording followed by two sessions on the track (a "value" and a "risk" session) and another rest recording.
-===== Exercise =====
+Let's focus on the first rest segment, using the ''ExpKeys'' to restrict the data:
-Recall the plots you made in [[analysis:nsb2015:week2|Module 2]]. The piece of data you zoomed in on contained sharp wave-ripple complexes (SWRs) recorded from the dorsal CA1 area of the hippocampus, manifest as brief (~200ms) high-frequency oscilation in the LFP  (the "ripple") which often rides on top of a slower deflection (the "sharp wave"), like this:
+<code matlab>
+% restrict to prerecord, leaving some time (10s) before rat actually goes on track
+csc_pre = restrict(csc,0,csc.cfg.ExpKeys.TimeOnTrack(1)-10);
+</code>
-{{ :analysis:course:hcsw.png?600 |}}
+Check if we have regular sampling with no gaps:
-These LFP events are associated with the synchronous activation of many cells, which is often structured to form "replay", the sequential activation of place cells corresponding to a coherent spatial trajectory.
+<code matlab>
+% check if sampling is ok
+plot(diff(csc_preR)); % only minimal differences
+Fs = 1./mean(diff(csc_preR));
+</code>
-Studies of replay start with the detection of potential replay events. To do this, we need to isolate those features of the LFP that are associated with SWR events and distinguish them from those resulting from artifacts associated with chewing, grooming, et cetera.
+The slightly odd structure of the diffs is because of Neuralynx's 512-sample-per-block format. So, technically we don't have uniformly spaced samples, but it's close enough that we don't have to bother with interpolating.
-Based on the filtering concepts above, we can implement a workflow for detecting SWR events, as follows:
+Let's decimate to speed up subsequent processing:
 <code matlab>
-%% some hippocampus data
+dsf = 4;
-cd('C:\data\R042-2013-08-18_recording');
+csc_pre.data = decimate(csc_pre.data,dsf);
-cfg = []; cfg.fc = {'R042-2013-08-18-CSC03a.ncs'};
+csc_pre.tvec = downsample(csc_pre.tvec,dsf);
-lfp = LoadCSC(cfg);
+csc_pre.cfg.hdr{1}.SamplingFrequency = csc_pre.cfg.hdr{1}.SamplingFrequency./dsf;
+</code>
-%% filter in SWR band
+(Note that the above is a quick hack of something that should really be implemented in a nice function, because the above commands really correspond to one conceptual operation of decimation, presenting an opportunity for making the code more readable and maintainable. But I haven't gotten around to that yet.)
-cfg = [];
-cfg.f = [140 220];
-cfg.display_filter = 0;
-SWRf = FilterLFP(cfg,lfp);
+Now we can compute the spectrum:
-%% obtain power and z-score it
+<code matlab>
-SWRp = LFPpower([],SWRf);
+wSize = 1024;
-SWRp_z = zscore_tsd(SWRp);
+[Pxx,F] = periodogram(csc_pre.data,hamming(length(csc_pre.data)),length(csc_pre.data),Fs);
+plot(F,10*log10(Pxx),'k'); xlabel('Frequency (Hz)'); ylabel('Power (dB)');
+xlim([0 150]);
+</code>
-%% detect events
+This gives:
-cfg = [];
-cfg.method = 'raw';
-cfg.threshold = 3;
-cfg.dcn =  '>'; % return intervals where threshold is exceeded
-cfg.merge_thr = 0.05; % merge events closer than this
-cfg.minlen = 0.05; % minimum interval length
-SWR_evt = TSDtoIV(cfg,SWRp_z);
+{{ :analysis:course:week4_fig7.png?600 |}}
-%% to each event, add a field with the max z-scored power (for later selection)
+Note that we are plotting not the raw power spectrum but rather the ''10*log10()'' of it; this is conventional so that the resulting units are in decibels (dB), the unit of signal power also applied to sound waves.
-cfg = [];
-cfg.method = 'max'; % 'min', 'mean'
-cfg.label = 'maxSWRp'; % what to call this in iv, i.e. usr.label
-SWR_evt = AddTSDtoIV(cfg,SWR_evt,SWRp_z);
+Regardless, it doesn't look very good! The estimates look very noisy, a characteristic of the periodogram method.
-%% select only those events of >5 z-scored power
+☛ Edit the above code to compute a Welch spectrum instead, using a Hamming window of the same size and 50% overlap. Much better!
-cfg = [];
-cfg.dcn = '>';
-cfg.threshold = 5;
-SWR_evt = SelectIV(cfg,SWR_evt);
+This PSD exhibits a number of features characteristic for neural data:
+  - The sharp spike at 60Hz is //line noise//, an artifact of the 60Hz oscillation in the North American alternating-current power supply.
+  - The overall shape of the PSD is described as ''1/f'', i.e. to a first approximation, the signal power is proportional to the inverse of the frequency. This is characteristic of many biological signals, see e.g. the Buzsaki (2006) book for a discussion of the possible significance of this.
+  - This LFP is recorded from the ventral striatum, which has been reported to exhibit delta (~4Hz), theta (~7-9Hz), beta (~12-20Hz), and gamma (~40-100Hz) oscillations. These are are all apparent as "humps" protruding from the underlying 1/f shape.
-%% plot events in highlighted on top of full lfp
+==== Exercises ====
-PlotTSDfromIV([],SWR_evt,lfp);
-%% ..or the events alone (fixed 200ms window centered at event time)
+Enter your code in a ''sandbox.m'' file using Cell Mode as usual.
-close all;
-cfg = [];
+☛ 1. Compute the PSD of "white noise", i.e. a random signal where each sample is drawn independently from the open interval (0,1) with equal probability. Is it 1/f? How would you describe its shape? //Hint//: use the MATLAB function ''rand()''.
-cfg.display = 'iv';
-cfg.mode = 'center';
-cfg.fgcol = 'k';
-PlotTSDfromIV(cfg,SWR_evt,lfp);
+☛ 2. Compute the PSD of a LFP, simultaneously recorded with the signal above but now from the hippocampus. The ''ExpKeys'' specify the filename of a verified hippocampal LFP (in the ''GoodSWR'' field, for "good sharp wave-ripple complexes", characteristic of the hippocampus). How does it compare to the ventral striatal LFP?
-%% ..hold on (highlight edges of event on top of previous plot)
-cfg = [];
-cfg.display = 'iv';
-cfg.fgcol = 'r';
-PlotTSDfromIV(cfg,SWR_evt,lfp);
+☛ 3. For both LFPs, explore the effects of the window size parameter of the Welch power spectrum, making sure to vary it across an order of magnitude. What do you notice?
-</code>
-☛ Try it, and inspect the results. What strategies can you think of to evaluate the accuracy and precision of the above detection? How might the workflow be improved?
+☛ 4. Compare the PSD following the use of ''decimate()'' as in the above example, to a PSD obtained from downsampling without using ''decimate()''. Are there any differences?

User Tools

Site Tools

Differences

Page Tools