Home > matlabmk > pwr_spec.m

# pwr_spec

## PURPOSE pwr_spec-Computers the power spectrum of a real valued time series using

## DESCRIPTION ``` 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:

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,
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}
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/2010```

## CROSS-REFERENCE INFORMATION This function calls:
This function is called by:

## SOURCE CODE ```0001 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:
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,
0059 end
0060
0061 if nargin<5,
0062     verblevel=2;
0063 end
0064
0065 [n_pts, n_trials]=size(data);
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```

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005