dftfilt2() - discrete complex wavelet filters
Usage:
>> wavelet = dftfilt2( freqs, cycles, srate, cyclefact)
Inputs:
freqs - frequency array
cycles - cycles array. If one value is given, 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
interpolation between these values for intermediate
frequencies
srate - sampling rate (in Hz)
cycleinc - ['linear'|'log'] increase mode if [min max] cycles is
provided in 'cycle' parameter. {default: 'linear'}
type - ['sinus'|'morlet'] wavelet type is a sinusoid with
cosine (real) and sine (imaginary) parts tapered by
a Hanning or Morlet function. 'morlet' is a typical Morlet
wavelet (with p=2*pi and sigma=0.7) best matching the
'sinus' Hanning taper) {default: 'morlet'}
Output:
wavelet - cell array of wavelet filters
Note: The length of the window is automatically computed from the
number of cycles and is always made odd.
Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/20030001 % dftfilt2() - discrete complex wavelet filters 0002 % 0003 % Usage: 0004 % >> wavelet = dftfilt2( freqs, cycles, srate, cyclefact) 0005 % 0006 % Inputs: 0007 % freqs - frequency array 0008 % cycles - cycles array. If one value is given, all wavelets have 0009 % the same number of cycles. If two values are given, the 0010 % two values are used for the number of cycles at the lowest 0011 % frequency and at the highest frequency, with linear 0012 % interpolation between these values for intermediate 0013 % frequencies 0014 % srate - sampling rate (in Hz) 0015 % 0016 % cycleinc - ['linear'|'log'] increase mode if [min max] cycles is 0017 % provided in 'cycle' parameter. {default: 'linear'} 0018 % type - ['sinus'|'morlet'] wavelet type is a sinusoid with 0019 % cosine (real) and sine (imaginary) parts tapered by 0020 % a Hanning or Morlet function. 'morlet' is a typical Morlet 0021 % wavelet (with p=2*pi and sigma=0.7) best matching the 0022 % 'sinus' Hanning taper) {default: 'morlet'} 0023 % Output: 0024 % wavelet - cell array of wavelet filters 0025 % 0026 % Note: The length of the window is automatically computed from the 0027 % number of cycles and is always made odd. 0028 % 0029 % Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003 0030 0031 % Copyright (C) 3/28/2003 Arnaud Delorme 8, SCCN/INC/UCSD, arno@salk.edu 0032 % 0033 % This program is free software; you can redistribute it and/or modify 0034 % it under the terms of the GNU General Public License as published by 0035 % the Free Software Foundation; either version 2 of the License, or 0036 % (at your option) any later version. 0037 % 0038 % This program is distributed in the hope that it will be useful, 0039 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0040 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0041 % GNU General Public License for more details. 0042 % 0043 % You should have received a copy of the GNU General Public License 0044 % along with this program; if not, write to the Free Software 0045 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 0046 0047 function wavelet = dftfilt2( freqs, cycles, srate, cycleinc, type); 0048 0049 if nargin < 3 0050 error('3 arguments required'); 0051 end; 0052 if nargin < 5 0053 type = 'morlet'; 0054 end; 0055 0056 % compute number of cycles at each frequency 0057 % ------------------------------------------ 0058 if length(cycles) == 1 0059 cycles = cycles*ones(size(freqs)); 0060 elseif length(cycles) == 2 0061 if nargin == 4 & strcmpi(cycleinc, 'log') % cycleinc 0062 cycles = linspace(log(cycles(1)), log(cycles(2)), length(freqs)); 0063 cycles = exp(cycles); 0064 else 0065 cycles = linspace(cycles(1), cycles(2), length(freqs)); 0066 end; 0067 end; 0068 0069 % compute wavelet 0070 for index = 1:length(freqs) 0071 0072 % number of cycles depend on window size 0073 % number of cycles automatically reduced if smaller window 0074 % note: as the number of cycle changes, the frequency shifts a little 0075 % this has to be fixed 0076 0077 winlen = cycles(index)*srate/freqs(index); 0078 winlenint = floor(winlen); 0079 if mod(winlenint,2) == 1, winlenint = winlenint+1; end; 0080 winval = linspace(winlenint/2, -winlenint/2, winlenint+1); 0081 0082 if strcmpi(type, 'sinus') % Hanning 0083 win = exp(2i*pi*freqs(index)*winval/srate); 0084 wavelet{index} = win .* hanning(length(winval))'; 0085 0086 else % Morlet 0087 t = freqs(index)*winval/srate; 0088 p = 2*pi; 0089 s = cycles(index)/5; 0090 wavelet{index} = exp(j*t*p)/sqrt(2*pi) .* ... 0091 (exp(-t.^2/(2*s^2))-sqrt(2)*exp(-t.^2/(s^2)-p^2*s^2/4)); 0092 end; 0093 end; 0094 0095 0096 return; 0097 0098 % testing 0099 % ------- 0100 wav1 = dftfilt2(5, 5, 256); wav1 = wav1{1}; 0101 abs1 = linspace(-floor(length(wav1)),floor(length(wav1)), length(wav1)); 0102 figure; plot(abs1, real(wav1), 'b'); 0103 0104 wav2 = dftfilt2(5, 3, 256); wav2 = wav2{1}; 0105 abs2 = linspace(-floor(length(wav2)),floor(length(wav2)), length(wav2)); 0106 hold on; plot(abs2, real(wav2), 'r'); 0107 0108 wav3 = dftfilt2(5, 1.4895990, 256); wav3 = wav3{1}; 0109 abs3 = linspace(-floor(length(wav3)),floor(length(wav3)), length(wav3)); 0110 hold on; plot(abs3, real(wav3), 'g'); 0111 0112 wav4 = dftfilt2(5, 8.73, 256); wav4 = wav4{1}; 0113 abs4 = linspace(-floor(length(wav4)),floor(length(wav4)), length(wav4)); 0114 hold on; plot(abs4, real(wav4), 'm'); 0115 0116 % more testing 0117 % ------------ 0118 freqs = exp(linspace(0,log(10),10)); 0119 win = dftfilt2(freqs, [3 3], 256, 'linear', 'sinus'); 0120 win = dftfilt2(freqs, [3 3], 256, 'linear', 'morlet'); 0121 0122 freqs = [12.0008 13.2675 14.5341 15.8007 17.0674 18.3340 19.6007 20.8673 22.1339 23.4006 24.6672 25.9339 27.2005 28.4671 29.7338 31.0004 32.2670 33.5337 34.8003 36.0670 37.3336 38.6002 39.8669 41.1335 42.4002 43.6668 44.9334 46.2001 47.4667 ... 0123 48.7334 50.0000]; 0124 0125 win = dftfilt2(freqs, [3 12], 256, 'linear'); size(win) 0126 0127 winsize = 0; 0128 for index = 1:length(win) 0129 winsize = max(winsize,length(win{index})); 0130 end; 0131 allwav = zeros(winsize, length(win)); 0132 for index = 1:length(win) 0133 wav1 = win{index}; 0134 abs1 = linspace(-(length(wav1)-1)/2,(length(wav1)-1)/2, length(wav1)); 0135 allwav(abs1+(winsize+1)/2,index) = wav1(:); 0136 end; 0137 figure; imagesc(imag(allwav)); 0138 0139 0140 % syemtric hanning function 0141 function w = hanning(n) 0142 if ~rem(n,2) 0143 w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1))); 0144 w = [w; w(end:-1:1)]; 0145 else 0146 w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1))); 0147 w = [w; w(end-1:-1:1)]; 0148 end