Nothing
## 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 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
## 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
kaiserord <- function(f, m, dev, Fs = 2) {
## parameter checking
if (length(f) != 2*length(m)-2)
stop("kaiserord must have one magnitude for each frequency band")
if (length(m) > 2 && any(m[1:(length(m)-2)] != m[3:length(m)]))
stop("kaiserord pass and stop bands must be strictly alternating")
if (length(dev) != length(m) && length(dev) != 1)
stop("kaiserord must have one deviation for each frequency band")
dev = min(dev)
if (dev <= 0)
stop("kaiserord must have dev > 0")
## use midpoints of the transition region for band edges
w = (f[seq(1, length(f), by = 2)] + f[seq(2, length(f), by = 2)])/Fs
## determine ftype
if (length(w) == 1)
if (m[1] > m[2]) ftype='low' else ftype='high'
else if (length(w) == 2)
if (m[1] > m[2]) ftype='stop' else ftype='pass'
else
if (m[1] > m[2]) ftype='DC-1' else ftype='DC-0'
## compute beta from dev
A = -20*log10(dev)
if (A > 50)
beta = 0.1102*(A-8.7)
else if (A >= 21)
beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21)
else
beta = 0.0
## compute n from beta and dev
dw = 2*pi*min(f[seq(2, length(f), by = 2)] - f[seq(1, length(f), by = 2)])/Fs
n = max(1, ceiling((A-8) / (2.285*dw)))
## if last band is high, make sure the order of the filter is even.
if ((m[1] > m[2]) == (length(w) %% 2 == 0) &&
n %% 2 == 1)
n = n+1
FilterOfOrder(n = n, Wc = w, type = ftype, beta = beta)
}
#!demo
###Fs = 11025
###op = par(mfrow=c(2,2), mar=c(3,3,1,1))
###for ( i in 1 : 4) {
### if (i==1) {
### bands=c(1200, 1500); mag=c(1, 0); dev=c(0.1, 0.1)
### } else if (i==2) {
### bands=c(1000, 1500); mag=c(0, 1); dev=c(0.1, 0.1)
### } else if (i==3) {
### bands=c(1000, 1200, 3000, 3500); mag=c(0, 1, 0); dev=0.1
### } else if (i==4) {
### bands=100*c(10, 13, 15, 20, 30, 33, 35, 40)
### mag=c(1, 0, 1, 0, 1); dev=0.05
### }
### kaisprm = kaiserord(bands, mag, dev, Fs)
### with(kaisprm, {
### d <<- max(1,trunc(n/10))
### if (mag[length(mag)]==1 && d %% 2 ==1)
### d<<-d+1
### f1 <<- freqz(fir1(n,w,ftype,kaiser(n+1,beta),'noscale'), Fs = Fs, plot = FALSE)
### f2 <<- freqz(fir1(n-d,w,ftype,kaiser(n-d+1,beta),'noscale'), Fs = Fs, plot = FALSE)
### })
### plot(f1$f,abs(f1$h), col = "blue", type="l", xlab = "", ylab = "")
### lines(f2$f,abs(f2$h), col = "red")
### legend("right", paste("order", c(kaisprm$n-d, kaisprm$n)), col=c("red", "blue"), lty=1, bty="n")
### b = c(0, bands, Fs/2)
### for ( i in seq(2,length(b),by=2)) {
### hi=mag[i/2]+dev[1]; lo=max(mag[i/2]-dev[1],0)
### lines(c(b[i-1], b[i], b[i], b[i-1], b[i-1]),c(hi, hi, lo, lo, hi))
### }
###}
###par(op)
#! %--------------------------------------------------------------
#! % 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 end demo to show detail at criteria box corners
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.