## Copyright (C) 2018 Ionescu Vlad ## ## 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 3 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, see . ## -*- texinfo -*- ## @deftypefn {Function File} address@hidden =} fir (@var{n}, @var{f}) ## @deftypefnx {Function File} address@hidden =} fir (@var{n}, @var{f}, @var{a}) ## @deftypefnx {Function File} address@hidden =} fir (@var{n}, @var{f}, @var{a}, @var{win}) ## @deftypefnx {Function File} address@hidden =} fir (@var{n}, @var{f}, @var{a}, @var{win}, @var{scale}) ## ## Create a single, or multiband FIR filter of order @var{n} with the frequency ## band edges @var{f}, thus @var{h} will have a length of @var{n}+1. For types ## I & II, specifying either a band edge at DC, or one at Nyquist, but not both, ## will mean the difference between having a DC, or not. For types III & IV, the ## same difference will be by either specifying both DC and Nyquist, or none. ## ## @var{a} specifies the gain of each passband. If @var{f} starts at 0, the ## first passband is between 0 and f(1), the 2nd between f(2) and f(3), if they ## exist, and so on. Similar for @var{f} with 1 at the end. ## ## An optional shaping @var{win} can be given as a vector with length @var{n}+1. ## Defaults to a rectangular (boxcar) window of length @var{n}+1. ## ## @var{scale} is a flag with the value "scale", or "y", meaning the filter's ## coefficients will be normalized such that the magnitude response at DC, or ## the Nyquist, or the center of the first passband is 1, depending on the case. ## ## To apply the filter, use the return vector @var{h} with the @code{filter} ## function, for example @code{y = filter (h, 1, x)}. ## ## Examples: ## @example ## freqz (fir (40, [0 0.37])); ## freqz (fir (15, [0.41 1], "y")); ## ... ## @end example ## @seealso{filter, fir1, fir2} ## @end deftypefn function h = fir (N, F, varargin) ## First check the required number of arguments if (nargin < 2 || nargin > 5) print_usage; endif ## Order must be a one one element vector... if (length (N) != 1) error ("The order (N) must be a vector of size 1.") endif ## ...and proper valued if ((N <= 0) || ischar (N)) error ("The order (N) must be positive definite.") endif ## The frequencies must be positive values, ... if (min (F) < 0) error ("The frequencies must be positive.") endif ## ... strictly increasing values, ... if (min (diff (F) <= 0)) error ("The frequencies must be strictly increasing.") endif ## ... and their range between 0 and 1, closed interval at both ends. if ((F(1) < 0) || (F(end) > 1)) error ("The frequencies must lie in the interval [0..1], with 1 being Nyquist.") endif ## Only real numbers if (iscomplex (F)) error ("The values of the frequency vector must not be imaginary.") endif ## Default values N = fix (N); # silently consider the integer part of N lenF = length(F); fType = (F(1) == 0 && F(end) == 1) || (F(1) > 0 && F(end) < 1); scale = 0; win = ones (1, N+1); A = ones(1, ceil (lenF/2)); ## Out of the 3 possible variable arguments, 1 will be string for k=1:length (varargin) if (isempty (varargin)) continue; # as per fir1.m, line 78, this may be to circumvent a bug endif switch (varargin{k}) case {"scale" "y"} scale = 1; otherwise if (ischar (varargin{k})) error ("The flag must be \"scale\", or \"y\".") endif if (length (varargin{k}) == N+1 ) win = varargin{k}; else A = varargin{k}; endif endswitch endfor ## Some helpers DC = (F(1) == 0); lenA = length (A); oddF = mod(lenF, 2); endF = lenF - (DC == 0 && oddN); oddN = mod (N+1, 2); midN = ceil ((N + 1)/2); M = N/2; W = F*pi; ## Handle errors from varargin ## A needs to have different lengths, depending on the type of FIR if (! fType) # types I & II if (DC && lenA != ceil (lenF/2)) error ("For type I & II with DC, the length of the amplitude vector must be ceil(length(F)/2).") elseif (! DC && lenA != floor (lenF/2)) error ("For type I & II without DC, the length of the amplitude vector must be floor(length(F)/2).") endif else if (DC && lenA != lenF/2) error ("For type III & IV, the length of the amplitude vector must be length(F)/2.") endif endif ## Only real numbers if (iscomplex (A)) error ("The values of the amplitude vector must not be imaginary.") endif ## window if (length (win) != N+1) error ("The length of the window must be the same as the filter's.") endif ## Make sure the vectors are columns F = F(:)'; ## F needs some adjusting if ( (DC && ! fType && oddF) || (! DC && fType && oddF) ) F = [F 1]; lenF += 1; oddF = 1 - oddF; endif win = win(:)'; A = [A(:)' ; A(:)'](:)'; ## A, too if (! DC && ! fType && oddF) A = [A 0]; endif ## Begin carnage n = -M:1:M; h = zeros (1, N+1); ## Types I & II if (! fType) ## sin(0) or sin(pi) are useless, just account for the sign for k=1+DC:endF h += (-1)^k*A(k)*sin(W(k)*n)./(pi*n).*win; endfor ## Even N will result in h(midN)=inf. Here, all the frequencies are needed. if (oddN) h(midN) = 0; for k=1:lenF h(midN) += (-1)^k*A(k)*F(k); endfor endif ## Types III & IV else ## Types III & IV need to use (1+cos(x))/x for DC, cos(x)/x for no DC for k=1:lenF h += (-1)^(k + 1)*A(k)*cos(W(k)*n)./(pi*n).*win; endfor if (oddN) h(midN) = 0; # midpoint is null for type III endif endif ## Normalizing for unity gain if (scale) ## Not valid for types III & IV, one, because there can be no DC, and two, ## their sum is 0, so they are considered, by all practical means, bandpass ## and/or bandstop variants, which means the normalizing is done at the ## center of the first bandwidth, whether F(1)=0 or not. The same is valid ## for type II with Nyquist (so to speak). if (DC && ! fType) h = A(1)*h/sum (h); elseif (! DC && ! fType && oddN) h = A(end)*h/abs (polyval (h, -1)); else h = A(1)*h/abs (polyval (h, exp(-0.5i*(W(1) + W(2))) )); endif endif endfunction %% expected output %!test %! x = [0.006248855546642409; ... %! -0.0151891541716538; ... %! -0.02000749470712845; ... %! 3.703551396872747*10^-4; ... %! 0.02352639910729612; ... %! 0.01987710435709456; ... %! -0.0107156702183056; ... %! -0.03346910762203006; ... %! -0.01647493965701502; ... %! 0.02756343124616879; ... %! 0.0468059024239567; ... %! 0.006350821104335874; ... %! -0.06144311077699433; ... %! -0.07272786646455694; ... %! 0.02972318687964276; ... %! 0.2090466916579446; ... %! 0.3495187814185787; ... %! 0.3495187814185787; ... %! 0.2090466916579446; ... %! 0.02972318687964276; ... %! -0.07272786646455694; ... %! -0.06144311077699433; ... %! 0.006350821104335874; ... %! 0.0468059024239567; ... %! 0.02756343124616879; ... %! -0.01647493965701502; ... %! -0.03346910762203006; ... %! -0.0107156702183056; ... %! 0.01987710435709456; ... %! 0.02352639910729612; ... %! 3.703551396872747*10^-4; ... %! -0.02000749470712845; ... %! -0.0151891541716538; ... %! 0.006248855546642409]; %! N = 33; %! f = [0 0.37]; %! h = fir (N, f)'; %! assert (x, h, 1e-15); %% expected output %!test %! x = [-0.006511935135150087; ... %! -0.02855131730142876; ... %! -0.009709917741985558; ... %! 0.02887609738977625; ... %! 0.03026378946243989; ... %! -0.0178273752838507; ... %! -0.05195697062909502; ... %! -0.009830953088030569; ... %! 0.07107858475501896; ... %! 0.0692657030840666; ... %! -0.08418362416801213; ... %! -0.3017430389202977; ... %! 0.5824183786740378; ... %! -0.3017430389202977; ... %! -0.08418362416801213; ... %! 0.0692657030840666; ... %! 0.07107858475501896; ... %! -0.009830953088030569; ... %! -0.05195697062909502; ... %! -0.0178273752838507; ... %! 0.03026378946243989; ... %! 0.02887609738977625; ... %! -0.009709917741985558; ... %! -0.02855131730142876; ... %! -0.006511935135150087]; %! N = 24; %! f = [0.41 1]; %! h = fir (N, f, 'y')'; %! assert (x, h, 1e-15);