analysis:course-w16:week5
====== Differences ====== This shows you the differences between two versions of the page.
| Both sides previous revision Previous revision Next revision | Previous revision | ||
|
analysis:course-w16:week5 [2016/01/21 20:01] mvdm [Decomposing and reconstructing a signal] |
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]] | ||
| 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 215: | 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 258: | 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 302: | 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 348: | 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 480: | 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 496: | 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 510: | 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 521: | 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.1453406490.txt.gz · Last modified: 2018/04/17 15:21 (external edit)
Page Tools
