# v_zoomfft

## PURPOSE

V_ZOOMFFT DTFT evaluated over a linear frequency range Y=(X,N,M,S,D)

## SYNOPSIS

function [y,f]=v_zoomfft(x,n,m,s,d)

## DESCRIPTION

    V_ZOOMFFT    DTFT evaluated over a linear frequency range Y=(X,N,M,S,D)
     Inputs:
        x    vector (or matrix)
        n    reciprocal of normalized frequency increment (can be non-integer).
             The frequency increment is fs/n where fs is the sample frequency
             [default n=size(x,d)]
        m    mumber of output points is floor(m) [default m=n]
        s    starting frequency index (can be non-integer).
             The starting frequency is s*fs/n [default s=0]
        d    dimension along which to do fft [default d=first non-singleton]

     Outputs:
        y       Output dtft coefficients. y has the same dimensions as x except
                that size(y,d)=floor(m).
        f(1,m)  normalized frequencies (1 corresponds to fs)

     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 DTFT 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 v_zoomfft(x,n,n,0,d) is equivalent to fft(x,n,d) for n>=length(x).

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

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