User Tools

Site Tools


analysis:course-w16:week5

Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
analysis:course-w16:week5 [2016/01/21 14:37]
mvdm
analysis:course-w16:week5 [2018/07/07 10:19] (current)
Line 1: Line 1:
 ~~DISCUSSION~~ ~~DISCUSSION~~
- 
-:!: **UNDER CONSTRUCTION,​ DO NOT USE!** :!: 
  
 ===== Fourier series, transforms, power spectra ===== ===== Fourier series, transforms, power spectra =====
Line 14: Line 12:
 Resources: Resources:
  
-  * (bakcground) {{:​analysis:​course:​weeks_ch4_sinusoids.pdf|Chapter 4}} from Weeks, "​Digital Signal Processing Using MATLAB & Wavelets"​+  * (background) {{:​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"​   * (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]]   * (optional) [[http://​www.mathworks.com/​help/​signal/​ug/​spectral-analysis.html|MATLAB'​s guide to spectral analysis]]
Line 20: Line 18:
 ==== Introductory remarks ==== ==== Introductory remarks ====
  
-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.+From very early electrical recordings of the human brain (e.g. 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.
  
 ==== Generating and plotting basic sinusoids ==== ==== Generating and plotting basic sinusoids ====
Line 26: Line 24:
 To study oscillations we require //​periodic//​ functions of time, such as ''​sin(t)'',​ which repeat its values in regular intervals. To study oscillations we require //​periodic//​ functions of time, such as ''​sin(t)'',​ which repeat its values in regular intervals.
  
-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:+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:
  
 <code matlab> <code matlab>
Line 35: Line 33:
  
 f = 2; % frequency of sine to plot f = 2; % frequency of sine to plot
-y = sin(2*pi*f*tvec);​ % note sin() expects arguments in radians, not degrees (see sind())+y = sin(2*pi*f*tvec);​ % note sin() expects arguments in radians, not degrees (see also ''​sind()''​)
  
 stem(tvec,​y);​ stem(tvec,​y);​
Line 43: Line 41:
  
 <code matlab> <code matlab>
- 
 phi = pi/2; phi = pi/2;
- +  
-figure;+subplot(221)
 y = sin(2*pi*f*tvec + phi); % a phase shift y = sin(2*pi*f*tvec + phi); % a phase shift
 stem(tvec,​y);​ stem(tvec,​y);​
 hold on; hold on;
-plot(tvec,​cos(2*pi*f*tvec),'​LineWidth',​2) % notice, cosine is simply phase shifted sine +plot(tvec,​cos(2*pi*f*tvec),'​r--'​,'​LineWidth',​2) % notice, cosine is simply phase shifted sine 
-legend('​Sin (Phase Shifted)', 'Cos')+legend('​sin (phase-shifted)', 'cos'); 
-hold off+ 
 a = 2; a = 2;
- +subplot(222)
-figure;+
 y = a.*sin(2*pi*f*tvec + phi); % amplitude change y = a.*sin(2*pi*f*tvec + phi); % amplitude change
-stem(tvec,​y);​ +stem(tvec,​y); ​% note scale of y axis!
 </​code>​ </​code>​
  
Line 82: Line 76:
 </​code>​ </​code>​
  
-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.+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, referred to "​cross-frequency coupling":​ in this example, the phase of one (slower, "​message"​) signal is related ("​coupled"​) to the frequency of another (faster, "​carrier"​) signal.
  
 ==== Sums of sinusoids and harmonic series ==== ==== Sums of sinusoids and harmonic series ====
Line 95: Line 89:
  
 <code matlab> <code matlab>
-%%Harmonic ​series example+%% harmonic ​series example
 mag = [0.1 0 1.3 0.5]; % magnitudes for each term mag = [0.1 0 1.3 0.5]; % magnitudes for each term
 pha = [-pi/6 0 pi 2*pi/3]; % phases for each term pha = [-pi/6 0 pi 2*pi/3]; % phases for each term
 f = 2; % base frequency f = 2; % base frequency
 + 
 signal_out = zeros(size(tvec));​ signal_out = zeros(size(tvec));​
 for ii = 1:​numel(mag) % note, the book chapter uses i, not best practice! for ii = 1:​numel(mag) % note, the book chapter uses i, not best practice!
 + 
     this_signal = mag(ii)*cos(2*pi*f*ii*tvec + pha(ii));     this_signal = mag(ii)*cos(2*pi*f*ii*tvec + pha(ii));
     plot(tvec,​this_signal,'​r:'​);​ hold on;     plot(tvec,​this_signal,'​r:'​);​ hold on;
     signal_out = signal_out + this_signal;​ % build the sum     signal_out = signal_out + this_signal;​ % build the sum
-    ​+ 
 end end
-figure; 
 plot(tvec,​signal_out,'​LineWidth',​2);​ plot(tvec,​signal_out,'​LineWidth',​2);​
 </​code>​ </​code>​
Line 114: Line 107:
 ☛ 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''​...) ☛ 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''​...)
  
-It looks like we can create some interesting signals by summing sinusoids.+It looks like we can create some interesting signals ​(blue) ​by summing ​simple ​sinusoids ​(thin red lines)!
  
 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. 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.
Line 123: Line 116:
  
 <code matlab> <code matlab>
 +%%
 +rng('​default'​);​ % reset random number generator to reproducible state, so your plot will look like mine!
 +
 x = round(rand(1,​8)*10);​ % generate a length 8 vector of integers between 0 and 10 x = round(rand(1,​8)*10);​ % generate a length 8 vector of integers between 0 and 10
 xlen = length(x); xlen = length(x);
 + 
 % get magnitudes and phases of Fourier series % get magnitudes and phases of Fourier series
 X = fft(x); X = fft(x);
 Xmag = abs(X); % magnitudes, a_n Xmag = abs(X); % magnitudes, a_n
 Xphase = angle(X); % phases, phi_n Xphase = angle(X); % phases, phi_n
 + 
 n = 0:xlen-1; n = 0:xlen-1;
 t = 0:​0.05:​xlen-1;​ % a finer timescale to show the smooth signal later t = 0:​0.05:​xlen-1;​ % a finer timescale to show the smooth signal later
 + 
 for iH = xlen-1:-1:0 % reconstruct each harmonic 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;​     s(iH+1,:) = Xmag(iH+1)*cos(2*pi*n*iH/​xlen + Xphase(iH+1))/​xlen;​
Line 139: Line 135:
     % detail: xlen appears here because the fundamental frequency used by fft() depends on this     % detail: xlen appears here because the fundamental frequency used by fft() depends on this
 end end
- +  
-ssum = sum(s); +ssum = sum(s); ​% coarse timescale (original points) 
-smsum = sum(sm); +smsum = sum(sm); ​% fine timescale (to see full signal) 
 + 
 figure; figure;
 plot(n, x, '​go',​ t, smsum, '​b',​ n, ssum, '​r*'​);​ plot(n, x, '​go',​ t, smsum, '​b',​ n, ssum, '​r*'​);​
Line 148: Line 144:
 </​code>​ </​code>​
  
-This gives something like:+You should get:
  
-{{ :​analysis:​course:​week4_fig1.png?nocache&600 }}+{{ :​analysis:​course-w16:fourier.png?nolink&​600 ​|}}
  
 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. 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.
Line 156: Line 152:
 ☛ 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.) ☛ 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.)
  
-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. +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 square of the magnitude spectrum ​is referred to as the signal //power// at different frequencies.
 ==== Interpreting the output of MATLAB'​s fft() function ====  ==== Interpreting the output of MATLAB'​s fft() function ==== 
  
Line 218: Line 213:
 This means that our signal is now no longer an integer number of periods. The resulting two-sided spectrum is: This means that our signal is now no longer an integer number of periods. The resulting two-sided spectrum is:
  
-{{ :​analysis:​course:​week4_fig3.png?600 |}}+{{ :​analysis:​course-w16:fourier2.png?nolink&600 |}}
  
 Notice that: Notice that:
Line 261: Line 256:
 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. 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.
  
-FIXME Demonstrate that you get the same effect by zero-padding (rather than making the signal longer)+☛ Demonstrate that you get the same effect by zero-padding ​the original 21-point signal ​(rather than making the signal longer, as in the above example).
  
 ==== Spectral leakage ==== ==== Spectral leakage ====
Line 305: Line 300:
 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. ​ 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. ​
  
-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.+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. ​(If the idea of convolution is new to you, you can think of it as a kind of "​blurring":​ for instance, convolving a signal with a 5-point rectangular pulse is equivalent to taking the running average with a 5-point moving window. For a graphical exploration of this idea, see [[http://​pages.jh.edu/​~signals/​convolve/​ | this website]] or look at the MATLAB doc for the ''​conv()''​ function.) ​
  
 Thus, it becomes important to understand the spectral properties of the windows used for Fourier analysis. Let's plot a few: Thus, it becomes important to understand the spectral properties of the windows used for Fourier analysis. Let's plot a few:
Line 351: Line 346:
 ☛ 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. ☛ 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.
  
-FIXME Clarify if the integral of the windows should be 1 (to preserve power estimates) and where this gets done in the code.+Note that the integral of the windows should be 1to preserve power estimates, as guaranteed by the normalization statement ​in the code above.
 ==== Robust spectral estimation methods ===== ==== Robust spectral estimation methods =====
  
Line 483: Line 478:
 <code matlab> <code matlab>
 % check if sampling is ok % check if sampling is ok
-plot(diff(csc_preR)); % only minimal differences +plot(diff(csc_pre.tvec)); % only minimal differences 
-Fs = 1./​mean(diff(csc_preR));+Fs = 1./​mean(diff(csc_pre.tvec));
 </​code>​ </​code>​
  
Line 499: Line 494:
  
 (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.) (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.)
 +
 +☛ Now that you have decimated the data, also update your sampling frequency ''​Fs''​.
  
 Now we can compute the spectrum: Now we can compute the spectrum:
Line 513: Line 510:
 {{ :​analysis:​course:​week4_fig7.png?​600 |}} {{ :​analysis:​course:​week4_fig7.png?​600 |}}
  
-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. ​+Note that we are plotting not the raw power spectrum but rather the ''​10*log10()''​ of it; this is a convention similar to that of the definition of the decibel ​(dB), the unit of signal power also applied to sound waves. The values on your vertical axis may look different, depending on what version of ''​LoadCSC()''​ was used.
  
 Regardless, it doesn'​t look very good! The estimates look very noisy, a characteristic of the periodogram method. Regardless, it doesn'​t look very good! The estimates look very noisy, a characteristic of the periodogram method.
Line 524: Line 521:
   - 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.   - 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.
  
-==== Exercises ====+☛ 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?
  
-Enter your code in a ''​sandbox.m'' ​file using Cell Mode as usual.+☛ Compute the PSD of "white noise",​ i.e. random signal where each sample is drawn independently from the interval (0,1) with equal probability. Is it 1/f? How would you describe its shape? //Hint//: use the MATLAB function ​''​rand()''​.
  
-☛ 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 probabilityIs it 1/f? How would you describe its shape//Hint//: use the MATLAB function ''​rand()''​.+☛ For both LFPs, explore ​the effects ​of the window size parameter of the Welch power spectrummaking sure to vary it across an order of magnitudeWhat do you notice?
  
-☛ 2. Compute ​the PSD of a LFP, simultaneously recorded with the signal above but now from the hippocampus. The ''​ExpKeys'' ​specify ​the filename of 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?+☛ Compare ​the PSD following the use of ''​decimate()'' ​as in the above example, to PSD obtained from downsampling without using ''​decimate()''​. ​Are there any differences?
  
-☛ 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?+==== Challenge ====
  
-☛ 4. Compare the PSD following the use of ''​decimate()''​ as in the above exampleto a PSD obtained from downsampling without using ''​decimate()''​Are there any differences?+★ Compute PSDs for some piece of data you collected. Make sure you specify briefly ​in the code commentsor short separate document, what this data describes, and its basic properties such as sampling rate and any filtering that may have been applied. When computing the PSD, systematically vary a few of the parameters ​(zero-padding,​ window size, window type). Do you feel some are more suitable than othersWhy? Also note some of your conclusions or interpretations based on what the PSDs show.
analysis/course-w16/week5.1453405028.txt.gz · Last modified: 2018/07/07 10:19 (external edit)