Mercurial > hg > octave-nkf

comparison scripts/signal/fftfilt.m @ 1026:9fc405c8c06c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-01-11 21:17:01 by jwe]
author jwe
date Wed, 11 Jan 1995 21:17:01 +0000
parents 3470f1e25a79
children 611d403c7f3d
comparison
equal deleted inserted replaced
1025:f558749713f1 1026:9fc405c8c06c
1 # Copyright (C) 1995 John W. Eaton
2 #
3 # This file is part of Octave.
4 #
5 # Octave is free software; you can redistribute it and/or modify it
6 # under the terms of the GNU General Public License as published by the
7 # Free Software Foundation; either version 2, or (at your option) any
8 # later version.
9 #
10 # Octave 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
13 # for more details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with Octave; see the file COPYING. If not, write to the Free
17 # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
18
1 function y = fftfilt (b, x, N) 19 function y = fftfilt (b, x, N)
2 20
3 # usage: fftfilt (b, x [, N]) 21 # usage: fftfilt (b, x [, N])
4 # 22 #
5 # y = fftfilt (b, x) filters x with the FIR filter b using the FFT. 23 # y = fftfilt (b, x) filters x with the FIR filter b using the FFT.
6 # y = fftfilt (b, x, N) uses the overlap-add method to filter x with 24 # y = fftfilt (b, x, N) uses the overlap-add method to filter x with
7 # b using an N-point FFT. 25 # b using an N-point FFT.
8 26
9 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994 27 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
10 # Copyright Dept of Statistics and Probability Theory TU Wien
11 28
12 # Reference: Oppenheim & Schafer (1989). Discrete-time Signal 29 # Reference: Oppenheim & Schafer (1989). Discrete-time Signal
13 # Processing (Chapter 8). Prentice-Hall. 30 # Processing (Chapter 8). Prentice-Hall.
14 31
15 # If N is not specified explicitly, we do not use the overlap-add 32 # If N is not specified explicitly, we do not use the overlap-add
22 usage (" fftfilt (b, x [, N])"); 39 usage (" fftfilt (b, x [, N])");
23 endif 40 endif
24 41
25 [r_x, c_x] = size (x); 42 [r_x, c_x] = size (x);
26 [r_b, c_b] = size (b); 43 [r_b, c_b] = size (b);
27 if !( (min ([r_x c_x]) == 1) || (min ([r_b c_b]) == 1) ) 44 if (! (min ([r_x c_x]) == 1 || min ([r_b c_b]) == 1))
28 error ("fftfilt: both x and b should be vectors."); 45 error ("fftfilt: both x and b should be vectors");
29 endif 46 endif
30 l_x = r_x * c_x; 47 l_x = r_x * c_x;
31 l_b = r_b * c_b; 48 l_b = r_b * c_b;
32 49
33 if ((l_x == 1) && (l_b == 1)) 50 if ((l_x == 1) && (l_b == 1))
34 y = b * x; 51 y = b * x;
35 return; 52 return;
36 endif 53 endif
37 54
38 x = reshape (x, 1, l_x); 55 x = reshape (x, 1, l_x);
39 b = reshape (b, 1, l_b); 56 b = reshape (b, 1, l_b);
40 57
41 if (nargin == 2) 58 if (nargin == 2)
42 # use FFT with the smallest power of 2 which is >= length (x) + 59 # Use FFT with the smallest power of 2 which is >= length (x) +
43 # length (b) - 1 as number of points ... 60 # length (b) - 1 as number of points ...
44 N = 2^(ceil (log (l_x + l_b - 1) / log(2))); 61 N = 2^(ceil (log (l_x + l_b - 1) / log(2)));
45 y = ifft (fft (x, N) .* fft(b, N)); 62 y = ifft (fft (x, N) .* fft(b, N));
46 else 63 else
47 # use overlap-add method ... 64 # Use overlap-add method ...
48 if !(is_scalar (N)) 65 if !(is_scalar (N))
49 error ("fftfilt: N has to be a scalar"); 66 error ("fftfilt: N has to be a scalar");
50 endif 67 endif
51 N = 2^(ceil (log (max ([N l_b])) / log(2))); 68 N = 2^(ceil (log (max ([N l_b])) / log(2)));
52 L = N - l_b + 1; 69 L = N - l_b + 1;
53 B = fft (b, N); 70 B = fft (b, N);
54 R = ceil (l_x / L); 71 R = ceil (l_x / L);
55 y = zeros (1, l_x); 72 y = zeros (1, l_x);
56 for r=1:R; 73 for r = 1:R;
57 lo = (r - 1) * L + 1; 74 lo = (r - 1) * L + 1;
58 hi = min (r * L, l_x); 75 hi = min (r * L, l_x);
59 tmp = ifft (fft (x(lo:hi), N) .* B); 76 tmp = ifft (fft (x(lo:hi), N) .* B);
60 hi = min (lo+N-1, l_x); 77 hi = min (lo+N-1, l_x);
61 y(lo:hi) = y(lo:hi) + tmp(1:(hi-lo+1)); 78 y(lo:hi) = y(lo:hi) + tmp(1:(hi-lo+1));
62 endfor 79 endfor
63 endif 80 endif
64 81
65 y = reshape (y(1:l_x), r_x, c_x); 82 y = reshape (y(1:l_x), r_x, c_x);
66 83
67 # final cleanups: if both x and b are real respectively integer, y 84 # Final cleanups: if both x and b are real respectively integer, y
68 # should also be so 85 # should also be
69 if !(any (imag (x)) || any (imag (b))) 86
87 if (! (any (imag (x)) || any (imag (b))))
70 y = real (y); 88 y = real (y);
71 endif 89 endif
72 if !(any (x - round (x)) || any (b - round (b))) 90 if (! (any (x - round (x)) || any (b - round (b))))
73 y = round (y); 91 y = round (y);
74 endif 92 endif
75 93
76 endfunction 94 endfunction
77 95