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).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 General Public License as published by 0038 % the Free Software Foundation; either version 2 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 General Public License for more details. 0045 % 0046 % You can obtain a copy of the GNU General Public License from 0047 % http://www.gnu.org/copyleft/gpl.html or by writing to 0048 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 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;