Permalink
Fetching contributors…
Cannot retrieve contributors at this time
109 lines (98 sloc) 3.14 KB
%firws() - Designs windowed sinc type I linear phase FIR filter
%
% Usage:
% >> b = firws(m, f);
% >> b = firws(m, f, w);
% >> b = firws(m, f, t);
% >> b = firws(m, f, t, w);
%
% Inputs:
% m - filter order (mandatory even)
% f - vector or scalar of cutoff frequency/ies (-6 dB;
% pi rad / sample)
%
% Optional inputs:
% w - vector of length m + 1 defining window {default hamming}
% t - 'high' for highpass, 'stop' for bandstop filter {default low-/
% bandpass}
%
% Output:
% b - filter coefficients
%
% Example:
% fs = 500; cutoff = 0.5; df = 1;
% m = firwsord('hamming', fs, df);
% b = firws(m, cutoff / (fs / 2), 'high', windows('hamming', m + 1));
%
% References:
% Smith, S. W. (1999). The scientist and engineer's guide to digital
% signal processing (2nd ed.). San Diego, CA: California Technical
% Publishing.
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
% firwsord, invfirwsord, kaiserbeta, windows
%123456789012345678901234567890123456789012345678901234567890123456789012
% Copyright (C) 2005 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%
% $Id$
function [b, a] = firws(m, f, t, w)
a = 1;
if nargin < 2
error('Not enough input arguments');
end
if length(m) > 1 || ~isnumeric(m) || ~isreal(m) || mod(m, 2) ~= 0 || m < 2
error('Filter order must be a real, even, positive integer.');
end
f = f / 2;
if any(f <= 0) || any(f >= 0.5)
error('Frequencies must fall in range between 0 and 1.');
end
if nargin < 3 || isempty(t)
t = '';
end
if nargin < 4 || isempty(w)
if ~isempty(t) && ~ischar(t)
w = t;
t = '';
else
w = windows('hamming', (m + 1));
end
end
w = w(:)'; % Make window row vector
b = fkernel(m, f(1), w);
if length(f) == 1 && strcmpi(t, 'high')
b = fspecinv(b);
end
if length(f) == 2
b = b + fspecinv(fkernel(m, f(2), w));
if isempty(t) || ~strcmpi(t, 'stop')
b = fspecinv(b);
end
end
% Compute filter kernel
function b = fkernel(m, f, w)
m = -m / 2 : m / 2;
b(m == 0) = 2 * pi * f; % No division by zero
b(m ~= 0) = sin(2 * pi * f * m(m ~= 0)) ./ m(m ~= 0); % Sinc
b = b .* w; % Window
b = b / sum(b); % Normalization to unity gain at DC
% Spectral inversion
function b = fspecinv(b)
b = -b;
b(1, (length(b) - 1) / 2 + 1) = b(1, (length(b) - 1) / 2 + 1) + 1;