dftfilt

 dftfilt() - discrete Fourier filter

 Usage:
   >> b = dftfilt(n,W,c,k,q)

 Inputs:
   n - number of input samples
   W - maximum angular freq. relative to n, 0 < W <= .5
   c - cycles
   k - oversampling
   q - [0;1] 0->fft, 1->c cycles

 Authors: Sigurd Enghoff, Arnaud Delorme & Scott Makeig, 
          SCCN/INC/UCSD, La Jolla, 8/1/98

0001 % dftfilt() - discrete Fourier filter
0002 %
0003 % Usage:
0004 %   >> b = dftfilt(n,W,c,k,q)
0005 %
0006 % Inputs:
0007 %   n - number of input samples
0008 %   W - maximum angular freq. relative to n, 0 < W <= .5
0009 %   c - cycles
0010 %   k - oversampling
0011 %   q - [0;1] 0->fft, 1->c cycles
0012 %
0013 % Authors: Sigurd Enghoff, Arnaud Delorme & Scott Makeig,
0014 %          SCCN/INC/UCSD, La Jolla, 8/1/98
0015 
0016 % Copyright (C) 8/1/98 Sigurd Enghoff & Scott Makei, SCCN/INC/UCSD, scott@sccn.ucsd.edu
0017 %
0018 % This program is free software; you can redistribute it and/or modify
0019 % it under the terms of the GNU General Public License as published by
0020 % the Free Software Foundation; either version 2 of the License, or
0021 % (at your option) any later version.
0022 %
0023 % This program is distributed in the hope that it will be useful,
0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0026 % GNU General Public License for more details.
0027 %
0028 % You should have received a copy of the GNU General Public License
0029 % along with this program; if not, write to the Free Software
0030 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0031 
0032 % 01-25-02 reformated help & license -ad
0033 
0034 % future developments
0035 % -------------------
0036 % input into dftfilt:
0037 % - lowfreq and maxfreq (of interest)
0038 % - lowcycle and maxcyle (ex: 3 cycles at low freq and 10 cycles at maxfreq)
0039 % - the delta in frequency: ex 0.5 Hz
0040 % The function should: compute the number of points (len) automatically
0041 % Warning with FFT compatibility
0042 % Still, something has to be done about the masking so that it would be comaptible
0043 
0044 function b = dftfilt(len,maxfreq,cycle,oversmp,wavfact)
0045 
0046 count = 1;
0047 for index = 1:1/oversmp:maxfreq*len/cycle % scan frequencies
0048     w(:,count) = j * index * cycle * linspace(-pi+2*pi/len, pi-2*pi/len, len)'; % exp(-w) is a sinus curve
0049     count = count+1; % -2*pi/len ensures that we really scan from -pi to pi without redundance (-pi=+pi)
0050 end;
0051 b = exp(-w);
0052 
0053 %srate = 2*pi/len;                            % Angular increment.
0054 %w = j * cycle * [0:srate:2*pi-srate/2]';    % Column.
0055 %x = 1:1/oversmp:maxfreq*len/cycle;            % Row.
0056 %b = exp(-w*x);                                % Exponentiation of outer product.
0057 
0058 for i = 1:size(b,2),
0059     m  = round(wavfact*len*(i-1)/(i+oversmp-1));    % Number of elements to discard.
0060     mu = round(m/2);                                % Number of upper elemnts.
0061     ml = m-round(m/2);                                % Number of lower elemnts.
0062     b(:,i) = b(:,i) .* [zeros(mu,1) ; hanning(len-m) ; zeros(ml,1)];
0063 end
0064 
0065 % syemtric hanning function
0066 function w = hanning(n)
0067 if ~rem(n,2)
0068    w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
0069    w = [w; w(end:-1:1)];
0070 else
0071    w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
0072    w = [w; w(end-1:-1:1)];
0073 end
0074