dftfilt3() - discrete complex wavelet filters
Usage:
>> [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin)
Inputs:
freqs - vector of frequencies of interest.
cycles - cycles array. If cycles=0, then the Hanning tapered Short-term FFT is used.
If one value is given and cycles>0, all wavelets have
the same number of cycles. If two values are given, the
two values are used for the number of cycles at the lowest
frequency and at the highest frequency, with linear or
log-linear interpolation between these values for intermediate
frequencies
srate - sampling rate (in Hz)
Optional Inputs: Input these as 'key/value pairs.
'cycleinc' - ['linear'|'log'] increase mode if [min max] cycles is
provided in 'cycle' parameter. {default: 'linear'}
'winsize' Use this option for Hanning tapered FFT or if you prefer to set the length of the
wavelets to be equal for all of them (e.g., to set the
length to 256 samples input: 'winsize',256). {default: [])
Note: the output 'wavelet' will be a matrix and it may be
incompatible with current versions of timefreq and newtimef.
'timesupport' The number of temporal standard deviation used for wavelet lengths {default: 7)
Output:
wavelet - cell array or matrix of wavelet filters
timeresol - temporal resolution of Morlet wavelets.
freqresol - frequency resolution of Morlet wavelets.
Note: The length of the window is always made odd.
Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003
Rey Ramirez, SCCN/INC/UCSD, La Jolla, 9/26/20060001 % dftfilt3() - discrete complex wavelet filters 0002 % 0003 % Usage: 0004 % >> [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin) 0005 % 0006 % Inputs: 0007 % freqs - vector of frequencies of interest. 0008 % cycles - cycles array. If cycles=0, then the Hanning tapered Short-term FFT is used. 0009 % If one value is given and cycles>0, all wavelets have 0010 % the same number of cycles. If two values are given, the 0011 % two values are used for the number of cycles at the lowest 0012 % frequency and at the highest frequency, with linear or 0013 % log-linear interpolation between these values for intermediate 0014 % frequencies 0015 % srate - sampling rate (in Hz) 0016 % 0017 % Optional Inputs: Input these as 'key/value pairs. 0018 % 'cycleinc' - ['linear'|'log'] increase mode if [min max] cycles is 0019 % provided in 'cycle' parameter. {default: 'linear'} 0020 % 'winsize' Use this option for Hanning tapered FFT or if you prefer to set the length of the 0021 % wavelets to be equal for all of them (e.g., to set the 0022 % length to 256 samples input: 'winsize',256). {default: []) 0023 % Note: the output 'wavelet' will be a matrix and it may be 0024 % incompatible with current versions of timefreq and newtimef. 0025 % 'timesupport' The number of temporal standard deviation used for wavelet lengths {default: 7) 0026 % 0027 % Output: 0028 % wavelet - cell array or matrix of wavelet filters 0029 % timeresol - temporal resolution of Morlet wavelets. 0030 % freqresol - frequency resolution of Morlet wavelets. 0031 % 0032 % Note: The length of the window is always made odd. 0033 % 0034 % Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003 0035 % Rey Ramirez, SCCN/INC/UCSD, La Jolla, 9/26/2006 0036 0037 % Copyright (C) 3/28/2003 Arnaud Delorme 8, SCCN/INC/UCSD, arno@salk.edu 0038 % 0039 % This program is free software; you can redistribute it and/or modify 0040 % it under the terms of the GNU General Public License as published by 0041 % the Free Software Foundation; either version 2 of the License, or 0042 % (at your option) any later version. 0043 % 0044 % This program is distributed in the hope that it will be useful, 0045 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0046 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0047 % GNU General Public License for more details. 0048 % 0049 % You should have received a copy of the GNU General Public License 0050 % along with this program; if not, write to the Free Software 0051 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 0052 0053 % 0054 % Revision 1.12 2006/09/25 rey r 0055 % Almost complete rewriting of dftfilt2.m, changing both Morlet and Hanning 0056 % DFT to be more in line with conventional implementations. 0057 % 0058 % Revision 1.11 2006/09/07 19:05:34 scott 0059 % further clarified the Morlet/Hanning distinction -sm 0060 % 0061 % Revision 1.10 2006/09/07 18:55:15 scott 0062 % clarified window types in help msg -sm 0063 % 0064 % Revision 1.9 2006/05/05 16:17:36 arno 0065 % implementing cycle array 0066 % 0067 % Revision 1.8 2004/03/04 19:31:03 arno 0068 % email 0069 % 0070 % Revision 1.7 2004/02/25 01:45:55 arno 0071 % sinus test 0072 % 0073 % Revision 1.6 2004/02/15 22:23:08 arno 0074 % implementing morlet wavelet 0075 % 0076 % Revision 1.5 2003/05/09 20:55:10 arno 0077 % adding hanning function 0078 % 0079 % Revision 1.4 2003/04/29 16:02:54 arno 0080 % header typos 0081 % 0082 % Revision 1.3 2003/04/29 01:09:16 arno 0083 % debug imaginary part 0084 % 0085 % Revision 1.2 2003/04/28 23:01:13 arno 0086 % *** empty log message *** 0087 % 0088 % Revision 1.1 2003/04/28 22:46:49 arno 0089 % Initial revision 0090 % 0091 0092 function [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin); 0093 0094 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0095 % Rey fixed all input parameter sorting. 0096 if nargin < 3 0097 error(' A minimum of 3 arguments is required'); 0098 end; 0099 numargin=length(varargin); 0100 if rem(numargin,2) 0101 error('There is an uneven number key/value inputs. You are probably missing a keyword or its value.') 0102 end 0103 varargin(1:2:end)=lower(varargin(1:2:end)); 0104 0105 % Setting default parameter values. 0106 cycleinc='linear'; 0107 winsize=[]; 0108 timesupport=7; % Setting default of 7 temporal standard deviations for wavelet's length. 0109 0110 for n=1:2:numargin 0111 keyword=varargin{n}; 0112 if strcmpi('cycleinc',keyword) 0113 cycleinc=varargin{n+1}; 0114 elseif strcmpi('winsize',keyword) 0115 winsize=varargin{n+1}; 0116 if ~mod(winsize,2) 0117 winsize=winsize+1; % Always set to odd length wavelets and hanning windows; 0118 end 0119 elseif strcmpi('timesupport',keyword) 0120 timesupport=varargin{n+1}; 0121 else 0122 error(['What is ' keyword '? The only legal keywords are: type, cycleinc, winsize, or timesupport.']) 0123 end 0124 end 0125 if isempty(winsize) & cycles==0 0126 error('If you are using a Hanning tapered FFT, please supply the winsize input-pair.') 0127 end 0128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0129 0130 0131 % compute number of cycles at each frequency 0132 % ------------------------------------------ 0133 type='morlet'; 0134 if length(cycles) == 1 & cycles(1)~=0 0135 cycles = cycles*ones(size(freqs)); 0136 elseif length(cycles) == 2 0137 if strcmpi(cycleinc, 'log') % cycleinc 0138 cycles = linspace(log(cycles(1)), log(cycles(2)), length(freqs)); 0139 cycles = exp(cycles); 0140 %cycles=logspace(log10(cycles(1)),log10(cycles(2)),length(freqs)); %rey 0141 else 0142 cycles = linspace(cycles(1), cycles(2), length(freqs)); 0143 end; 0144 end; 0145 if cycles==0 0146 type='sinus'; 0147 end 0148 0149 sp=1/srate; % Rey added this line (i.e., sampling period). 0150 % compute wavelet 0151 for index = 1:length(freqs) 0152 fk=freqs(index); 0153 if strcmpi(type, 'morlet') % Morlet. 0154 sigf=fk/cycles(index); % Computing time and frequency standard deviations, resolutions, and normalization constant. 0155 sigt=1./(2*pi*sigf); 0156 A=1./sqrt(sigt*sqrt(pi)); 0157 timeresol(index)=2*sigt; 0158 freqresol(index)=2*sigf; 0159 if isempty(winsize) % bases will be a cell array. 0160 tneg=[-sp:-sp:-sigt*timesupport/2]; 0161 tpos=[0:sp:sigt*timesupport/2]; 0162 t=[fliplr(tneg) tpos]; 0163 psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t)); 0164 wavelet{index}=psi; % These are the wavelets with variable number of samples based on temporal standard deviations (sigt). 0165 else % bases will be a matrix. 0166 tneg=[-sp:-sp:-sp*winsize/2]; 0167 tpos=[0:sp:sp*winsize/2]; 0168 t=[fliplr(tneg) tpos]; 0169 psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t)); 0170 wavelet(index,:)=psi; % These are the wavelets with the same length. 0171 % This is useful for doing time-frequency analysis as a matrix vector or matrix matrix multiplication. 0172 end 0173 elseif strcmpi(type, 'sinus') % Hanning 0174 tneg=[-sp:-sp:-sp*winsize/2]; 0175 tpos=[0:sp:sp*winsize/2]; 0176 t=[fliplr(tneg) tpos]; 0177 win = exp(2*i*pi*fk*t); 0178 wavelet(index,:) = win .* hanning(winsize)'; 0179 %wavelet{index} = win .* hanning(winsize)'; 0180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0181 end; 0182 end; 0183 0184 0185 0186 % symmetric hanning function 0187 function w = hanning(n) 0188 if ~rem(n,2) 0189 w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1))); 0190 w = [w; w(end:-1:1)]; 0191 else 0192 w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1))); 0193 w = [w; w(end-1:-1:1)]; 0194 end