zoomfft

ZOOMFFT    DFT evaluated over a linear frequency range Y=(X,N,M,S,D)
 Inputs:
    x    vector (or matrix)
    n    frequency increment is fs/n [default n=length(x)]
    m    mumber of output points is floor(m) [default m=n]
    s    starting frequency is s*fs/n [default s=0]
    d    dimension along which to do fft [default d=first non-singleton]

 This routine allows the evaluation of the DFT over an arbitrary range of
 frequencies; as its name implies this lets you zoom into a narrow portion
 of the spectrum.
 The DFT of X will be evaluated along dimension D at the M frequencies
 f=fs*(s+(0:m-1))/n where fs is the sample frequency. Note that N and S
 need not be integers although M will be rounded down to an integer.
 Thus zoomfft(x,n,n,0,d) is equivalent to fft(x,n,d) for n>=length(x).

0001 function y=zoomfft(x,n,m,s,d)
0002 %ZOOMFFT    DFT evaluated over a linear frequency range Y=(X,N,M,S,D)
0003 % Inputs:
0004 %    x    vector (or matrix)
0005 %    n    frequency increment is fs/n [default n=length(x)]
0006 %    m    mumber of output points is floor(m) [default m=n]
0007 %    s    starting frequency is s*fs/n [default s=0]
0008 %    d    dimension along which to do fft [default d=first non-singleton]
0009 %
0010 % This routine allows the evaluation of the DFT over an arbitrary range of
0011 % frequencies; as its name implies this lets you zoom into a narrow portion
0012 % of the spectrum.
0013 % The DFT of X will be evaluated along dimension D at the M frequencies
0014 % f=fs*(s+(0:m-1))/n where fs is the sample frequency. Note that N and S
0015 % need not be integers although M will be rounded down to an integer.
0016 % Thus zoomfft(x,n,n,0,d) is equivalent to fft(x,n,d) for n>=length(x).
0017 
0018 % [1] L.R.Rabiner,  R.W.Schafer and C.M.Rader, "The chirp z-transform algorithm"
0019 %     IEEE Trans. Audio Electroacoustics 17 (2), 86–92 (1969).
0020 
0021 %      Copyright (C) Mike Brookes 2007
0022 %      Version: $Id: zoomfft.m 713 2011-10-16 14:45:43Z dmb $
0023 %
0024 %   VOICEBOX is a MATLAB toolbox for speech processing.
0025 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0026 %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 %   This program is free software; you can redistribute it and/or modify
0029 %   it under the terms of the GNU General Public License as published by
0030 %   the Free Software Foundation; either version 2 of the License, or
0031 %   (at your option) any later version.
0032 %
0033 %   This program is distributed in the hope that it will be useful,
0034 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0035 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0036 %   GNU General Public License for more details.
0037 %
0038 %   You can obtain a copy of the GNU General Public License from
0039 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0040 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0042 persistent n0 k0 s0 m0 b c h g
0043 e=size(x);
0044 p=prod(e);
0045 if nargin<5
0046     d=find(e>1);
0047     if ~isempty(d)
0048         d=d(1);
0049     else
0050         d=1;
0051     end
0052 end
0053 k=e(d);
0054 q=p/k;
0055 if d==1
0056     z=reshape(x,k,q);
0057 else
0058     z=shiftdim(x,d-1);
0059     r=size(z);
0060     z=reshape(z,k,q);
0061 end
0062 if nargin<2 || isempty(n)
0063     n=k;
0064 end
0065 if nargin<3 || isempty(m)
0066     m=floor(n);
0067 else
0068     m=floor(m);
0069 end
0070 if nargin<4 || isempty(s)
0071     s=0;
0072 end
0073 l=pow2(nextpow2(m+k-1));    % round up to next power of 2
0074 if n==fix(n) && s==fix(s) && n<2*l && n>=k
0075     a=fft(z,n,1);           % quickest to do a normal fft
0076     y=a(1+mod(s:s+m-1,n),:);
0077 else
0078     % can precaluclate all this for fixed n, k, s and m
0079     if isempty(b) || n~=n0 || k~=k0 || s~=s0 || m~=m0
0080         n0=n;
0081         k0=k;
0082         s0=s;
0083         m0=m;
0084         b=exp(1i*pi*mod((s+(1-k:m-1)').^2,2*n)/n);
0085         c=conj(b(k:k+m-1));
0086         h=fft(b,l,1);
0087         g=exp(-1i*pi*mod(((0:k-1)').^2,2*n)/n);
0088     end
0089     a=ifft(fft(z.*repmat(g,1,q),l,1).*repmat(h,1,q)); % calculate correlation
0090     y=a(k:k+m-1,:).*repmat(c,1,q);
0091 end
0092 if d==1
0093     e(d)=m;
0094     y=reshape(y,e);
0095 else
0096     r(1)=m;
0097     y=shiftdim(reshape(y,r),length(e)+1-d);
0098 end