Diff of /fourier/chirpzt.m [971186] .. [b59441] Maximize Restore

  Switch to side-by-side view

--- a/fourier/chirpzt.m
+++ b/fourier/chirpzt.m
@@ -1,19 +1,85 @@
-function c = chirpzt(f,K,deltao,o,dim)
+function c = chirpzt(f,K,fdiff,foff,fs,dim)
 %CHIRPZT Chirped Z-transform
-%   Usage:  c = chirpzt(f,o0,deltao,K,dim)
+%   Usage:  c = chirpzt(f,K,fdiff)
+%           c = chirpzt(f,K,fdiff,foff)
+%           c = chirpzt(f,K,fdiff,foff,fs)
 %
 %   Input parameters:
 %         f      : Input data.
 %         K      : Number of values.
-%         deltao : Angle increment.
-%         o      : Starting angle. 
+%         fdiff  : Frequency increment.
+%         foff   : Starting frequency. 
+%         fs     : Sampling frequency. 
 %
 %   Output parameters:
 %         c      : Coefficient vector.
 %
-%   `c = chirpzt(f,o0,deltao,K)` computes *K* samples of the discrete-time 
-%   fourier transform DTFT *c* of *f* at values $c(k)=F(o0+k*deltao)$ 
-%   for $k=0,\dots,K-1$. The values are computed along dimension *dim*.
+%   `c = chirpzt(f,K,fdiff,foff)` computes *K* samples of the discrete-time 
+%   fourier transform DTFT *c* of *f* at values $c(k+1)=F(2*pi*(foff+k*fdiff))$
+%   for $k=0,\dots,K-1$ where $F=DTFT(f)$. Values `foff` and `fdiff` should
+%   be in range of $0-1$. If `foff` is ommited or empty, it is considered to
+%   be 0. If `fdiff` is ommited or empty, *K* equidistant values 
+%   $c(k+1)=F(2*pi*k/K)$ are computed. If even *K* is ommited or empty, 
+%   input length is used instead resulting in the same values as `fft` does.
+%
+%   `c = chirpzt(f,K,fdiff,foff,fs)` computes coefficients using frequency 
+%   values relative to *fs* $c(k+1)=F(2*pi*(foff+k*fdiff)/fs)$ for $k=0,\dots,K-1$.
+%
+%   The input *f* is processed along the first non-singleton dimension or
+%   along dimension *dim* if specified.
+%
+%   Examples:
+%   ---------
+%   
+%   Calculating DTFT samples of interest (aka zoom FFT):::
+% 
+%     % Generate input signal
+%     fs = 8000;
+%     L = 2^10;
+%     k = (0:L-1).';
+%     f1 = 400;
+%     f2 = 825;
+%     f = 5*sin(2*pi*k*f1/fs + pi/4) + 2*sin(2*pi*k*f2/fs - pi/3);
+%
+%     % This is equal to fft(f)
+%     ck = chirpzt(f,L);
+%
+%     %chirpzt to FFT error:
+%     norm(ck-fft(f))
+%
+%     % Frequency "resolution" in Hz
+%     fdiff = 0.4;
+%     % Frequency offset in Hz
+%     foff = 803.9;
+%     % Number of frequency values
+%     K = 125;
+%     % DTFT samples. The frequency range of interest is 803.9-853.5 Hz
+%     ckchzt = chirpzt(f,K,fdiff,foff,fs);
+%
+%     % Plot modulus of coefficients
+%     figure(1);
+%     fax=foff+fdiff.*(0:K-1);
+%     hold on;
+%     stem(k/L*fs,abs(ck),'k');
+%     stem(fax,abs(ckchzt),'r:');
+%     set(gca,'XLim',[foff,foff+K*fdiff]);
+%     set(gca,'YLim',[0 1065]);
+%     xlabel('f[Hz]');
+%     ylabel('|ck|');
+%
+%     % Plot phase of coefficients
+%     figure(2);
+%     hold on;
+%     stem(k/L*fs,angle(ck),'k');
+%     stem(fax,angle(ckchzt),'r:');
+%     set(gca,'XLim',[foff,foff+K*fdiff]);
+%     set(gca,'YLim',[-pi pi]);
+%     xlabel('f[Hz]');
+%     ylabel('angle(ck)');
+%
+%   See also: gga
+%
+%   References: raschra69
 
 %% Check the input arguments
 if nargin < 1
@@ -24,37 +90,46 @@
     error('%s: X must be a nonempty vector or a matrix.',upper(mfilename))
 end
 
-if nargin<5
+if nargin<6
   dim=[];  
 end;
 
 [f,~,Ls,~,dim,permutedsize,order]=assert_sigreshape_pre(f,[],dim,'CHIRPZT');
 
 if nargin > 1  && ~isempty(K)
-   if ~isreal(K) || ~isscalar(K)
-      error('%s: o0 must be a real scalar.',upper(mfilename))
+   if ~isreal(K) || ~isscalar(K) || rem(K,1)~=0
+      error('%s: K must be a real integer.',upper(mfilename))
    end
 else
    K = size(f,1);
 end
 
-if nargin > 2  && ~isempty(deltao)
+if nargin > 2  && ~isempty(fdiff)
    if ~isreal(K) || ~isscalar(K)
-      error('%s: deltao must be a real scalar.',upper(mfilename))
+      error('%s: fdiff must be a real scalar.',upper(mfilename))
    end
 else
-   deltao = 2*pi/K;
+   fdiff = 1/K;
 end
 
-if nargin > 3  && ~isempty(o)
+if nargin > 3  && ~isempty(foff)
    if ~isreal(K) || ~isscalar(K)
-      error('%s: o must be a real scalar.',upper(mfilename))
+      error('%s: foff must be a real scalar.',upper(mfilename))
    end
 else
-   o = 0;
+   foff = 0;
 end
 
-c = comp_chirpzt(f,K,deltao,o);
+if nargin > 4  && ~isempty(fs)
+   if ~isreal(fs) || ~isscalar(fs)
+      error('%s: fs must be a real scalar.',upper(mfilename))
+   end
+else
+   fs = 1;
+end
+
+
+c = comp_chirpzt(f,K,2*pi*fdiff/fs,2*pi*foff/fs);
 
 
 permutedsize(1)=K;