annotate scripts/signal/stft.m @ 6547:4fb053f24fd6

[project @ 2007-04-19 21:47:40 by jwe]
author jwe
date Thu, 19 Apr 2007 21:47:41 +0000
parents 34f96dd5441b
children 93c65f2a5668
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
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
7 ## the Free Software Foundation; either version 2, or (at your option)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
8 ## 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
3922
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
16 ## along with Octave; see the file COPYING. If not, write to the Free
5307
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 3922
diff changeset
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 3922
diff changeset
18 ## 02110-1301, USA.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
19
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
20 ## -*- texinfo -*-
6547
4fb053f24fd6 [project @ 2007-04-19 21:47:40 by jwe]
jwe
parents: 6046
diff changeset
21 ## @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
22 ## 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
23 ## @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
24 ## points and an increment of @var{inc} points.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
25 ##
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
26 ## 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
27 ## is applied:
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
28 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
29 ## @table @asis
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
30 ## @item hanning
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
31 ## w_type = 1
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
32 ## @item hamming
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
33 ## w_type = 2
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
34 ## @item rectangle
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
35 ## w_type = 3
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
36 ## @end table
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
37 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
38 ## 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
39 ##
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
40 ## 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
41 ## @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
42 ## @var{w_type} = 1.
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
43 ##
3449
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
44 ## @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
45 ## 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
46 ## frequencies.
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
47 ##
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
48 ## @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
49 ## 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
50 ## 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
51 ## synthesis function.
858695b3ed62 [project @ 2000-01-18 04:08:59 by jwe]
jwe
parents: 3426
diff changeset
52 ## @end deftypefn
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
53
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
54 ## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
55 ## Description: Short-term Fourier transform
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
56
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
57 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
58
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
59 ## default values of unspecified arguments
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
60 if (nargin < 5)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
61 w_type = 1;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
62 if (nargin < 4)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
63 coef = 64;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
64 if (nargin < 3)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
65 inc = 24;
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
66 if (nargin < 2)
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
67 win = 80;
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
68 endif
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
69 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
70 endif
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
71 elseif (nargin == 5)
5443
ec8c33dcd1bf [project @ 2005-09-08 01:40:57 by jwe]
jwe
parents: 5307
diff changeset
72 if (ischar (w_type))
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
73 if (strcmp (w_type, "hanning"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
74 w_type = 1;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
75 elseif (strcmp (w_type, "hamming"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
76 w_type = 2;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
77 elseif (strcmp (w_type, "rectangle"))
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
78 w_type = 3;
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
79 else
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
80 error ("stft: unknown window type `%s'", w_type);
3191
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 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
83 else
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5443
diff changeset
84 print_usage ();
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
85 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
86
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
87 ## check whether X is a vector
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
88 [nr, nc] = size (X);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
89 if (nc != 1)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
90 if (nr == 1)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
91 X = X';
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
92 nr = nc;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
93 else
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
94 error ("stft: X must be a vector");
3191
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 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
97
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
98 num_coef = 2 * coef;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
99 if (win > num_coef)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
100 win = num_coef;
3457
e031284eea27 [project @ 2000-01-19 08:49:56 by jwe]
jwe
parents: 3449
diff changeset
101 printf ("stft: window size adjusted to %f\n", win);
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
102 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
103 num_win = fix ((nr - win) / inc);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
104
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
105 ## compute the window coefficients
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
106 if (w_type == 3) # rectangular window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
107 WIN_COEF = ones (win, 1);
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
108 elseif (w_type == 2) # Hamming window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
109 WIN_COEF = hamming (win);
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
110 else # Hanning window
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
111 WIN_COEF = hanning (win);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
112 endif
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3191
diff changeset
113
3191
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
114 ## 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
115 Z = zeros (num_coef, num_win + 1);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
116 start = 1;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
117 for i = 0:num_win
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
118 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
119 start = start + inc;
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
120 endfor
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
121
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
122 Y = fft (Z);
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
123
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
124 if (nargout == 1)
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
125 Y = abs (Y(1:coef, :));
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
126 else
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
127 c = [win, inc, w_type];
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
128 endif
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
129
e4f4b2d26ee9 [project @ 1998-10-23 05:43:59 by jwe]
jwe
parents:
diff changeset
130 endfunction