X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fkaiserord.m;fp=octave_packages%2Fsignal-1.1.3%2Fkaiserord.m;h=2e7cf25fd3d9da23966f9b17b9ba0ae6221faf18;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/signal-1.1.3/kaiserord.m b/octave_packages/signal-1.1.3/kaiserord.m new file mode 100644 index 0000000..2e7cf25 --- /dev/null +++ b/octave_packages/signal-1.1.3/kaiserord.m @@ -0,0 +1,147 @@ +## Copyright (C) 2000 Paul Kienzle +## +## 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 . + +## usage: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs]) +## +## Returns the parameters needed for fir1 to produce a filter of the +## desired specification from a kaiser window: +## n: order of the filter (length of filter minus 1) +## Wn: band edges for use in fir1 +## beta: parameter for kaiser window of length n+1 +## ftype: choose between pass and stop bands +## b = fir1(n,Wn,kaiser(n+1,beta),ftype,'noscale'); +## +## f: frequency bands, given as pairs, with the first half of the +## first pair assumed to start at 0 and the last half of the last +## pair assumed to end at 1. It is important to separate the +## band edges, since narrow transition regions require large order +## filters. +## m: magnitude within each band. Should be non-zero for pass band +## and zero for stop band. All passbands must have the same +## magnitude, or you will get the error that pass and stop bands +## must be strictly alternating. +## dev: deviation within each band. Since all bands in the resulting +## filter have the same deviation, only the minimum deviation is +## used. In this version, a single scalar will work just as well. +## fs: sampling rate. Used to convert the frequency specification into +## the [0, 1], where 1 corresponds to the Nyquist frequency, fs/2. +## +## The Kaiser window parameters n and beta are computed from the +## relation between ripple (A=-20*log10(dev)) and transition width +## (dw in radians) discovered empirically by Kaiser: +## +## / 0.1102(A-8.7) A > 50 +## beta = | 0.5842(A-21)^0.4 + 0.07886(A-21) 21 <= A <= 50 +## \ 0.0 A < 21 +## +## n = (A-8)/(2.285 dw) +## +## Example +## [n, w, beta, ftype] = kaiserord([1000,1200], [1,0], [0.05,0.05], 11025); +## freqz(fir1(n,w,kaiser(n+1,beta),ftype,'noscale'),1,[],11025); + +## TODO: order is underestimated for the final test case: 2 stop bands. +## TODO: octave> ftest("kaiserord") # shows test cases + +function [n, w, beta, ftype] = kaiserord(f, m, dev, fs) + + if (nargin<2 || nargin>4) + print_usage; + endif + + ## default sampling rate parameter + if nargin<4, fs=2; endif + + ## parameter checking + if length(f)!=2*length(m)-2 + error("kaiserord must have one magnitude for each frequency band"); + endif + if any(m(1:length(m)-2)!=m(3:length(m))) + error("kaiserord pass and stop bands must be strictly alternating"); + endif + if length(dev)!=length(m) && length(dev)!=1 + error("kaiserord must have one deviation for each frequency band"); + endif + dev = min(dev); + if dev <= 0, error("kaiserord must have dev>0"); endif + + ## use midpoints of the transition region for band edges + w = (f(1:2:length(f))+f(2:2:length(f)))/fs; + + ## determine ftype + if length(w) == 1 + if m(1)>m(2), ftype='low'; else ftype='high'; endif + elseif length(w) == 2 + if m(1)>m(2), ftype='stop'; else ftype='pass'; endif + else + if m(1)>m(2), ftype='DC-1'; else ftype='DC-0'; endif + endif + + ## compute beta from dev + A = -20*log10(dev); + if (A > 50) + beta = 0.1102*(A-8.7); + elseif (A >= 21) + beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21); + else + beta = 0.0; + endif + + ## compute n from beta and dev + dw = 2*pi*min(f(2:2:length(f))-f(1:2:length(f)))/fs; + n = max(1,ceil((A-8)/(2.285*dw))); + + ## if last band is high, make sure the order of the filter is even. + if ((m(1)>m(2)) == (rem(length(w),2)==0)) && rem(n,2)==1, n = n+1; endif +endfunction + +%!demo +%! Fs = 11025; +%! for i=1:4 +%! if i==1, +%! subplot(221); bands=[1200, 1500]; mag=[1, 0]; dev=[0.1, 0.1]; +%! elseif i==2 +%! subplot(222); bands=[1000, 1500]; mag=[0, 1]; dev=[0.1, 0.1]; +%! elseif i==3 +%! subplot(223); bands=[1000, 1200, 3000, 3500]; mag=[0, 1, 0]; dev=0.1; +%! elseif i==4 +%! subplot(224); bands=100*[10, 13, 15, 20, 30, 33, 35, 40]; +%! mag=[1, 0, 1, 0, 1]; dev=0.05; +%! endif +%! [n, w, beta, ftype] = kaiserord(bands, mag, dev, Fs); +%! d=max(1,fix(n/10)); +%! if mag(length(mag))==1 && rem(d,2)==1, d=d+1; endif +%! [h, f] = freqz(fir1(n,w,ftype,kaiser(n+1,beta),'noscale'),1,[],Fs); +%! hm = freqz(fir1(n-d,w,ftype,kaiser(n-d+1,beta),'noscale'),1,[],Fs); +%! plot(f,abs(hm),sprintf("r;order %d;",n-d), ... +%! f,abs(h), sprintf("b;order %d;",n)); +%! b = [0, bands, Fs/2]; hold on; +%! for i=2:2:length(b), +%! hi=mag(i/2)+dev(1); lo=max(mag(i/2)-dev(1),0); +%! plot([b(i-1), b(i), b(i), b(i-1), b(i-1)],[hi, hi, lo, lo, hi],"c;;"); +%! endfor; hold off; +%! endfor +%! +%! %-------------------------------------------------------------- +%! % A filter meets the specifications if its frequency response +%! % passes through the ends of the criteria boxes, and fails if +%! % it passes through the top or the bottom. The criteria are +%! % met precisely if the frequency response only passes through +%! % the corners of the boxes. The blue line is the filter order +%! % returned by kaiserord, and the red line is some lower filter +%! % order. Confirm that the blue filter meets the criteria and +%! % the red line fails. + +%!# XXX FIXME XXX extend demo to show detail at criteria box corners