wiki

User Tools

Site Tools

Trace:
analysis:course-w16:week5

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

Link to this comparison view

Next revision 		Previous revision
analysis:course-w16:week5 [2016/01/21 19:36]
mvdm created 		analysis:course-w16:week5 [2018/04/17 15:20] (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]]
  
 ==== 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 1, to 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. a 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 probability. Is 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 spectrum, making sure to vary it across an order of magnitude. What 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 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?+☛ 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?
  
-☛ 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 example, to 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 comments, or a 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 others? Why? Also note some of your conclusions or interpretations based on what the PSDs show.
analysis/course-w16/week5.1453405006.txt.gz · Last modified: 2018/04/17 15:21 (external edit)

Page Tools

Donate Driven by DokuWiki Template From SyntaxSeed