Author: wmb
Date: 2009-08-04 03:32:33 +0200 (Tue, 04 Aug 2009)
New Revision: 1271
Modified:
cpu/x86/iirfilter.fth
Log:
Commented IIR upsampling filter code.
Modified: cpu/x86/iirfilter.fth
===================================================================
--- cpu/x86/iirfilter.fth 2009-08-04 00:46:03 UTC (rev 1270)
+++ cpu/x86/iirfilter.fth 2009-08-04 01:32:33 UTC (rev 1271)
@@ -1,12 +1,22 @@
-purpose: IIR digital filter support
+purpose: IIR digital upsampling filter
-\ iir-cascade implements the "Folded Cascade First Order Filters"
-\ structure as shown on page 274 of _Multirate Signal Processing
-\ for Communications Systems_ by fredric j harris, Prentice Hall 2004
-\ ISBN 0-13-146511-2
+\ This file implements x2 and x3 upsampling filters using techniques
+\ described in chapter 10 of _Multirate Signal Processing for Communications
+\ Systems_ by fredric j harris, Prentice Hall 2004, ISBN 0-13-146511-2
+\ The filters may be cascaded for x6 upsampling (e.g. for 8kHz -> 48kHz)
-d# 14 constant mulscale
+d# 15 constant mulscale \ Bits for fractional multipliers
+0 [if]
+
+\ iir-cascade-z2 implements a folded cascade of multiple sections
+\ of a second-order IIR filter block in a polynomial in Z^2.
+\ This is the basic building block for half-band all-pass
+\ recursive filters, as described in section 10.2 of the book
+\ cited above.
+\ This code is commented out because we don't use it directly,
+\ instead using the -z1 variant below in a polyphase structure.
+
code iir-cascade-z2 ( sample weights z #sections -- out )
cx pop \ CX: Number of sections
@@ -28,20 +38,29 @@
0 [bp] dx mov
dx 4 [bp] mov
- bx 0 [bp] mov \ Update M[i]
+ bx 0 [bp] mov \ Update Z[i]
- 8 [bp] bp lea \ Point to next M[i]
+ 8 [bp] bp lea \ Point to next Z[i]
loopa
0 [bp] dx mov
dx 4 [bp] mov
- ax 0 [bp] mov \ Update M[i]
+ ax 0 [bp] mov \ Update Z[i]
bp pop \ Restore BP
si pop \ Restore SI
ax push \ Return value
c;
+\ iir-cascade-z1 is the same basic structure as iir-cascade-z2,
+\ but the Z^-2 delays have been replaced by Z^-1. It is used
+\ in the polyphase upsampler "up2". When the rate is doubled,
+\ a zero sample is inserted between each pair of input samples,
+\ causing every other computation of the Z^-2 polynomial to be
+\ unnecessary. The basic block can thus be restructured to have
+\ a single delay element, and the block only has to run at half
+\ the rate. See section 10.6 (figure 10.48) of the book.
+
code iir-cascade-z1 ( sample weights z #sections -- out )
cx pop \ CX: Number of sections
@@ -54,25 +73,32 @@
bx push \ Saved BP on stack
begin
ax bx mov \ Save last value
- d# 4 [bp] ax sub \ last - M[i+1]
+ d# 4 [bp] ax sub \ last - Z[i+1]
- 0 [si] imul \ alpha[i] * (last - M[i+1]) (kills DX)
+ 0 [si] imul \ alpha[i] * (last - Z[i+1]) (kills DX)
4 [si] si lea \ Increment alpha pointer
mulscale # dx ax shrd \ Scale down by multiplier scale factor
- d# 0 [bp] ax add \ M[i] + alpha[i] * (last - M[i+1])
+ d# 0 [bp] ax add \ Z[i] + alpha[i] * (last - Z[i+1])
- bx 0 [bp] mov \ Update M[i]
+ bx 0 [bp] mov \ Update Z[i]
- 4 [bp] bp lea \ Point to next M[i]
+ 4 [bp] bp lea \ Point to next Z[i]
loopa
- ax 0 [bp] mov \ Update M[i]
+ ax 0 [bp] mov \ Update Z[i]
bp pop \ Restore BP
si pop \ Restore SI
ax push \ Return value
c;
+\ iir-cascade-z2g implements the generalized second-order
+\ low pass filter cascade as shown in section 10.5.1 (figure 10.35)
+\ of the book. The generalized form allows arbitrary bandwidths
+\ as opposed to the half-band restriction of the simplifed form.
+\ We use this form for 3x upsampling so we can set the lowpass
+\ frequency to Fnyquist/3.
+
code iir-cascade-z2g ( sample weights z #sections -- out )
cx pop \ CX: Number of sections
@@ -121,12 +147,14 @@
c;
+\ Convert a decimal fraction to a scaled integer multiplier
: mul: ( "coef" -- )
safe-parse-word push-decimal $number drop pop-base ( n )
1 mulscale lshift d# 1,000,000,000 */ ( n' )
,
;
+\ Clip too-large sample values to +- 15 bits
code saturate ( s -- s' )
ax pop
d# 32767 # ax cmp > if
@@ -146,44 +174,56 @@
;
[then]
-\ upsample by 4
-
2 constant #coefs
\ This is for a 9-pole, 9-zero 2-path all-pass halfband IIR filter
+\ From page 278 of the book
create weights0 mul: .101467517 mul: .612422841 \ Even path
create weights1 mul: .342095596 mul: .867647439 \ Odd path
0 [if]
\ This is for a 5-pole, 5-zero 2-path all-pass halfband IIR filter
+\ From page 276 of the book
create weights2 mul: .141348600 \ Even path
create weights3 mul: .589994800 \ Odd path
[then]
\ This is for a 5-pole, 5-zero 2-path all-pass third-band IIR filter
\ G0 and G1 use iir-cascade-z2g, H1 uses iir-cascade-z1
+\ This is derived from the 5-pole halfband coefficients according
+\ to the formulas on page 291, with fb/fs = 1/6 (fb/fnyquist = 1/3).
+\ The book has a missing sign somewhere, because without the minus
+\ signs on c1 and b, you end up with the corner frequency at 2/6
+\ instead of 1/6. I had to figure that out the hard way, by trial
+\ and error (looking at spectra of impulse responses). If you
+\ replace b with -b in equation 10.25, the c1's work out automatically.
+
+\ c1 c2
create coefg0 mul: -.605502011 mul: .211004022 \ Even path G0
create coefg1 mul: -.817448745 mul: .634897489 \ Odd path G1
+\ b
create coefh1 mul: -.267949192
-#coefs 1+ 2* /n* constant buflen
-buflen buffer: z0 \ First upsampler even path delay buffer
-buflen buffer: z1 \ First upsampler odd path delay buffer
-buflen buffer: z2 \ 3x upsampler even path G0 delay buffer
-buflen buffer: z3 \ 3x upsampler odd path G1 delay buffer
-buflen buffer: z4 \ 3x upsampler odd path H1 delay buffer
+2 1+ /n* constant buflen0 \ 2 stages of Z^-1 plus final feedback
-: init-upsample
- z0 buflen erase
- z1 buflen erase
- z2 buflen erase
- z3 buflen erase
- z4 buflen erase
+buflen0 buffer: z0 \ First upsampler even path delay buffer
+buflen0 buffer: z1 \ First upsampler odd path delay buffer
+
+1 1+ 2* /n* constant buflen0 \ 1 stage of Z^-2 plus final feedback
+
+buflen1 buffer: z2 \ 3x upsampler even path G0 delay buffer
+buflen1 buffer: z3 \ 3x upsampler odd path G1 delay buffer
+buflen1 buffer: z4 \ 3x upsampler odd path H1 delay buffer
+
+: init-upsample ( -- )
+ z0 buflen0 erase
+ z1 buflen0 erase
+ z2 buflen1 erase
+ z3 buflen1 erase
+ z4 buflen1 erase
;
-0 value lastsample
-
\ This is a polyphase decomposition of upsampling by 2
\ The filter prototype for each path is a polynomial in Z^2,
\ but when you apply the polyphase transform with Nobel's identity,
@@ -194,6 +234,8 @@
;
0 [if]
+0 value lastsample
+
\ This is the base rate (non-polyphase) version of up2
\ Both paths run for each sample, with Z^-1 between the
\ paths, summing the paths to get a lowpass output.
@@ -234,6 +276,7 @@
: 8khz>48khz ( inbuf #samples stride outbuf -- )
to sample-stride ( inbuf #samples outbuf )
sample-outp ! ( inbuf #samples )
+ init-upsample ( inbuf #samples )
0 ?do ( inbuf )
dup sample-stride + ( inbuf inbuf' )
swap <w@ ( inbuf' sample )
@@ -247,6 +290,7 @@
: 16khz>48khz ( inbuf #samples stride outbuf -- )
to sample-stride ( inbuf #samples outbuf )
sample-outp ! ( inbuf #samples )
+ init-upsample ( inbuf #samples )
0 ?do ( inbuf )
dup sample-stride + ( inbuf inbuf' )
swap <w@ ( inbuf' sample )
@@ -254,3 +298,28 @@
loop ( inbuf )
drop
;
+
+0 [if]
+\ Here is some Matlab/Octave code to do the "coefg.." coefficient
+\ transformations
+function retval = c1(a, b)
+ retval = 2 * b * (1 + a) / (1 + a * b * b);
+end
+function retval = c2(a, b)
+ retval = ((b * b) + a) / (1 + a * b * b);
+end
+function retval = cvtb(frac)
+ retval = (1 - tan(pi * frac)) / (1 + tan(pi * frac));
+end
+function cvtfilt (bw) % bw is fb/fs
+ a0 = .141348600
+ a1 = .589994800
+
+ b = cvtb(bw)
+ c10 = c1(a0, b)
+ c20 = c2(a0, b)
+ c11 = c1(a1, b)
+ c21 = c2(a1, b)
+end
+cvtfilt(1/6);
+[then]
