pwr_spec-Computers the power spectrum of a real valued time series using
one of various windows (i.e, tapers: boxcar, hann) and a fast fourier
transform. Note that power spectra are scaled such that if you use a
boxcar window, sum(pwr)=sum(data.^2). Other software for computing power
spectra may differ in how the scale the power spectra.
Usage:
>> [pwr, f]=pwr_spec(data,srate,wind,pad,verblevel);
Required Inputs:
data - [time x trial matrix] One or more trials of a time series
srate - [scalar] Sampling rate
Optional Inputs
wind - [string] Window function (i.e., "taper") applied to the time
series before computing the power spectra. Options are:
'boxcar': No taper {default}
'hann': Hann window
'hamming': Hamming window
pad - [integer] Padding factor for the FFT. Can take values -1, 0,
1, 2, etc... -1 corresponds to no padding. 0 corresponds to
no padding if the number of time points is a power of 2 or
padding to the next highest power of 2. 1 corresponds to
padding to the second next highest power of 2, etc... For
example, assuming there are 500 time points and if pad = -1,
we do not pad; if pad = 0, we pad the FFT to 512 points,
if pad=1, we pad to 1024 points etc... {default: 0}
verblevel - An integer specifiying the amount of information you want
this function to provide about what it is doing during runtime.
Options are:
0 - quiet, only show errors, warnings, and EEGLAB reports
1 - stuff anyone should probably know
2 - stuff you should know the first time you start working
with a data set {default value}
3 - stuff that might help you debug (show all
reports)
Outputs:
pwr - [frequency x trial matrix] Power at each frequency and trial.
First frequency is 0 and last frequency is the Nyquist frequency.
Power is scaled such that the sum(pwr)=sum(data.^2)
f - Frequencies (in Hz) corresponding to the first dimension of pwr
Author:
David Groppe
Kutaslab, 12/20100001 function [pwr, f]=pwr_spec(data,srate,wind,pad,verblevel) 0002 % pwr_spec-Computers the power spectrum of a real valued time series using 0003 % one of various windows (i.e, tapers: boxcar, hann) and a fast fourier 0004 % transform. Note that power spectra are scaled such that if you use a 0005 % boxcar window, sum(pwr)=sum(data.^2). Other software for computing power 0006 % spectra may differ in how the scale the power spectra. 0007 % 0008 % Usage: 0009 % >> [pwr, f]=pwr_spec(data,srate,wind,pad,verblevel); 0010 % 0011 % Required Inputs: 0012 % data - [time x trial matrix] One or more trials of a time series 0013 % srate - [scalar] Sampling rate 0014 % 0015 % Optional Inputs 0016 % wind - [string] Window function (i.e., "taper") applied to the time 0017 % series before computing the power spectra. Options are: 0018 % 'boxcar': No taper {default} 0019 % 'hann': Hann window 0020 % 'hamming': Hamming window 0021 % pad - [integer] Padding factor for the FFT. Can take values -1, 0, 0022 % 1, 2, etc... -1 corresponds to no padding. 0 corresponds to 0023 % no padding if the number of time points is a power of 2 or 0024 % padding to the next highest power of 2. 1 corresponds to 0025 % padding to the second next highest power of 2, etc... For 0026 % example, assuming there are 500 time points and if pad = -1, 0027 % we do not pad; if pad = 0, we pad the FFT to 512 points, 0028 % if pad=1, we pad to 1024 points etc... {default: 0} 0029 % verblevel - An integer specifiying the amount of information you want 0030 % this function to provide about what it is doing during runtime. 0031 % Options are: 0032 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0033 % 1 - stuff anyone should probably know 0034 % 2 - stuff you should know the first time you start working 0035 % with a data set {default value} 0036 % 3 - stuff that might help you debug (show all 0037 % reports) 0038 % 0039 % Outputs: 0040 % pwr - [frequency x trial matrix] Power at each frequency and trial. 0041 % First frequency is 0 and last frequency is the Nyquist frequency. 0042 % Power is scaled such that the sum(pwr)=sum(data.^2) 0043 % f - Frequencies (in Hz) corresponding to the first dimension of pwr 0044 % 0045 % Author: 0046 % David Groppe 0047 % Kutaslab, 12/2010 0048 0049 if nargin<2, 0050 error('You need at least two input arguments.'); 0051 end 0052 0053 if nargin<3, 0054 wind='boxcar'; 0055 end 0056 0057 if nargin<4, 0058 pad=0; 0059 end 0060 0061 if nargin<5, 0062 verblevel=2; 0063 end 0064 0065 [n_pts, n_trials]=size(data); 0066 if pad>=0, 0067 n_fft=2^(nextpow2(n_pts)+pad); 0068 else 0069 n_fft=n_pts; 0070 end 0071 0072 if verblevel>2, 0073 fprintf('Performing FFT on %d actual time points of data (%d after padding). There are %d trials.\n', ... 0074 n_pts,n_fft,n_trials); 0075 end 0076 0077 Nyq=srate/2; 0078 T=n_fft/srate; 0079 f=0:1/T:Nyq; 0080 n_freq=length(f); 0081 0082 pwr=zeros(n_freq,n_trials); %preallocate memory 0083 0084 %Compute vector to normalize FFT output such that sum of the squares of pwr 0085 %equals the sum of the time series squared. This involves dividing by 0086 %the number of time points and multiplying redundant frequencies by two 0087 %since we're only showing one member of redundant frequency pairs 0088 if rem(n_fft-1,2) 0089 nrmlize=[1 2*ones(1,n_freq-2) 1]'/n_pts; % 0 Hz and Nyquist are non redundant 0090 else 0091 nrmlize=[1 2*ones(1,n_freq-1)]'/n_pts; %only 0 Hz is non-redundant 0092 end 0093 0094 if isempty(wind) || strcmpi(wind,'boxcar'), 0095 if verblevel>1, 0096 fprintf('Using boxcar window (i.e., no window)\n'); 0097 end 0098 for tr=1:n_trials, 0099 Y=fft(data(:,tr),n_fft); 0100 %full_pwr=Y.*conj(Y)/n_pts; OLD 0101 full_pwr=Y.*conj(Y); 0102 pwr(:,tr)=full_pwr(1:n_freq).*nrmlize; 0103 end 0104 else 0105 if strcmpi(wind,'hann') || strcmpi(wind,'hanning'), 0106 if verblevel>1, 0107 fprintf('Using hann window\n'); 0108 end 0109 wind=hann(n_pts); 0110 elseif strcmpi(wind,'hamming'), 0111 if verblevel>1, 0112 fprintf('Using hamming window\n'); 0113 end 0114 wind=hamming(n_pts); 0115 end 0116 for tr=1:n_trials, 0117 Y=fft(data(:,tr).*wind,n_fft); 0118 %full_pwr=Y.*conj(Y)/n_pts; OLD 0119 full_pwr=Y.*conj(Y); 0120 pwr(:,tr)=full_pwr(1:n_freq).*nrmlize; 0121 end 0122 end