octave_packages/signal-1.1.3/welchwin.m

   1 ## Copyright (C) 2007 Muthiah Annamalai <muthiah.annamalai@uta.edu>
   2 ## Copyright (C) 2008-2009 Mike Gross <mike@appl-tech.com>
   3 ## Copyright (C) 2008-2009 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it under
   6 ## the terms of the GNU General Public License as published by the Free Software
   7 ## Foundation; either version 3 of the License, or (at your option) any later
   8 ## version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 ## details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License along with
  16 ## this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {[@var{w}] =} welchwin(@var{L},@var{c})
  20 ## Returns a row vector containing a Welch window, given by
  21 ## @var{w}(n)=1-(n/N-1)^2,   n=[0,1, ... @var{L}-1].
  22 ## Argument @var{L} is the length of the window.
  23 ## Optional argument @var{c} specifies a "symmetric" window (the default),
  24 ## or a "periodic" window.
  25 ##
  26 ## A symmetric window has zero at each end and maximum in the middle;
  27 ## @var{L} must be an integer larger than 2.
  28 ## @code{if c=="symmetric", N=(L-1)/2}
  29 ##
  30 ## A periodic window wraps around the cyclic interval [0,1, ... @var{L}-1],
  31 ## and is intended for use with the DFT  (functions fft(),
  32 ## periodogram() etc).
  33 ## @var{L} must be an integer larger than 1.
  34 ## @code{if c=="periodic", N=@var{L}/2}.
  35 ##
  36 ## @seealso{blackman, kaiser}
  37 ## @end deftypefn
  38
  39 function [w] = welchwin(L,c)
  40   if (nargin < 1 || nargin>2 )
  41     print_usage;
  42   endif
  43   symmetric=1;
  44   if ( nargin==2 && ! isempty(c) )
  45     if ( ! ischar(c) || size(c,1) != 1 ||
  46          ( ! strcmp(c,'periodic') && ! strcmp(c,'symmetric') ) )
  47       error( "arg 2 (c) must be \"periodic\" or \"symmetric\"" )
  48     endif
  49     symmetric = ! strcmp(c,'periodic');
  50   endif
  51   ##
  52   ## Periodic window is not properly defined if L<2.
  53   ## Symmetric window is not properly defined if L<3.
  54   min_L = 2 + symmetric;
  55   if ( ! isreal(L) || ! isscalar(L) || L<min_L || fix(L) != L )
  56     error("arg 1 (L) must be an integer larger than %d", min_L-1 );
  57   endif
  58   N = (L-symmetric)/2;
  59   n = 0:L-1;
  60   w = 1 - ((n-N)./N).^2;
  61 endfunction;
  62
  63 %!demo
  64 %! printf( "Short symmetric windows\n" );
  65 %! welchwin(5)
  66 %! welchwin(6)
  67 %! welchwin(6,"symmetric")
  68 %! welchwin(7)
  69 %! input( "Press ENTER to continue", "s" );
  70 %!
  71 %! printf( "Short periodic windows\n" );
  72 %! welchwin(6,"periodic")
  73 %! welchwin(7,"periodic")
  74 %! input( "Press ENTER to continue", "s" );
  75 %!
  76 %! L=32;
  77 %! printf( "Graph: single period of " );
  78 %! printf( "%d-point periodic (blue) and symmetric (red) windows\n", L );
  79 %! t=[0:L-1];
  80 %! xp=welchwin(L,"periodic");
  81 %! xs=welchwin(L,"symmetric");
  82 %! plot(t,xp,"b",t,xs,"r")
  83 %! input( "Press ENTER to continue", "s" );
  84 %!
  85 %! printf( "Graph: 4 periods of " );
  86 %! printf( "%d-point periodic (blue) and symmetric (red) windows\n", L );
  87 %! t2=[0:4*L-1];
  88 %! xp2=repmat(xp,1,4);
  89 %! xs2=repmat(xs,1,4);
  90 %! plot(t2,xp2,"b",t2,xs2,"r")
  91 %! input( "Press ENTER to continue", "s" );
  92 %!
  93 %! n=512;
  94 %! s=fftshift(max(1.0e-2,abs(fft(postpad(xp,n)))));
  95 %! f=[-0.5:1/n:0.5-1/n];
  96 %! printf( "%dx null-padded, power spectrum of %d-point window\n", n/L, L );
  97 %! semilogy(f,s)