annotate scripts/signal/stft.m @ 7016:93c65f2a5668

[project @ 2007-10-12 06:40:56 by jwe]
author jwe
date Fri, 12 Oct 2007 06:41:26 +0000
parents 4fb053f24fd6
children a1dbe9d80eee
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
1 ## Copyright (C) 1995, 1996, 1997 Andreas Weingessel
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
2 ##
3922
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
3 ## This file is part of Octave.
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
4 ##
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
6 ## under the terms of the GNU General Public License as published by
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
8 ## your option) any later version.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
9 ##
3922
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
13 ## General Public License for more details.
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
14 ##
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
17 ## <http://www.gnu.org/licenses/>.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
18
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
19 ## -*- texinfo -*-
6547
4fb053f24fd6 [project @ 2007-04-19 21:47:40 by jwe]
jwe
parents: 6046
diff changeset
20 ## @deftypefn {Function File} {[@var{y}, @var{c}] =} stft (@var{x}, @var{win_size}, @var{inc}, @var{num_coef}, @var{w_type})
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
21 ## Compute the short-term Fourier transform of the vector @var{x} with
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
22 ## @var{num_coef} coefficients by applying a window of @var{win_size} data
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
23 ## points and an increment of @var{inc} points.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
24 ##
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
25 ## Before computing the Fourier transform, one of the following windows
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
26 ## is applied:
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
27 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
28 ## @table @asis
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
29 ## @item hanning
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
30 ## w_type = 1
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
31 ## @item hamming
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
32 ## w_type = 2
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
33 ## @item rectangle
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
34 ## w_type = 3
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
35 ## @end table
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
36 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
37 ## The window names can be passed as strings or by the @var{w_type} number.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
38 ##
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
39 ## If not all arguments are specified, the following defaults are used:
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
40 ## @var{win_size} = 80, @var{inc} = 24, @var{num_coef} = 64, and
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
41 ## @var{w_type} = 1.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
42 ##
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
43 ## @code{@var{y} = stft (@var{x}, @dots{})} returns the absolute values
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
44 ## of the Fourier coefficients according to the @var{num_coef} positive
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
45 ## frequencies.
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
46 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
47 ## @code{[@var{y}, @var{c}] = stft (@code{x}, @dots{})} returns the
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
48 ## entire STFT-matrix @var{y} and a 3-element vector @var{c} containing
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
49 ## the window size, increment, and window type, which is needed by the
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
50 ## synthesis function.
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
51 ## @end deftypefn
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
52
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
53 ## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
54 ## Description: Short-term Fourier transform
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
55
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
56 function [Y, c] = stft(X, win, inc, coef, w_type)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
57
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
58 ## default values of unspecified arguments
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
59 if (nargin < 5)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
60 w_type = 1;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
61 if (nargin < 4)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
62 coef = 64;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
63 if (nargin < 3)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
64 inc = 24;
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
65 if (nargin < 2)
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
66 win = 80;
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
67 endif
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
68 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
69 endif
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
70 elseif (nargin == 5)
5443
ec8c33dcd1bf [project @ 2005-09-08 01:40:57 by jwe]
jwe
parents: 5307
diff changeset
71 if (ischar (w_type))
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
72 if (strcmp (w_type, "hanning"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
73 w_type = 1;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
74 elseif (strcmp (w_type, "hamming"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
75 w_type = 2;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
76 elseif (strcmp (w_type, "rectangle"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
77 w_type = 3;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
78 else
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
79 error ("stft: unknown window type `%s'", w_type);
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
80 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
81 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
82 else
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5443
diff changeset
83 print_usage ();
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
84 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
85
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
86 ## check whether X is a vector
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
87 [nr, nc] = size (X);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
88 if (nc != 1)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
89 if (nr == 1)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
90 X = X';
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
91 nr = nc;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
92 else
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
93 error ("stft: X must be a vector");
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
94 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
95 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
96
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
97 num_coef = 2 * coef;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
98 if (win > num_coef)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
99 win = num_coef;
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
100 printf ("stft: window size adjusted to %f\n", win);
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
101 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
102 num_win = fix ((nr - win) / inc);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
103
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
104 ## compute the window coefficients
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
105 if (w_type == 3) # rectangular window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
106 WIN_COEF = ones (win, 1);
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
107 elseif (w_type == 2) # Hamming window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
108 WIN_COEF = hamming (win);
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
109 else # Hanning window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
110 WIN_COEF = hanning (win);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
111 endif
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
112
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
113 ## create a matrix Z whose columns contain the windowed time-slices
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
114 Z = zeros (num_coef, num_win + 1);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
115 start = 1;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
116 for i = 0:num_win
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
117 Z(1:win, i+1) = X(start:start+win-1) .* WIN_COEF;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
118 start = start + inc;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
119 endfor
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
120
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
121 Y = fft (Z);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
122
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
123 if (nargout == 1)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
124 Y = abs (Y(1:coef, :));
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
125 else
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
126 c = [win, inc, w_type];
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
127 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
128
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
129 endfunction