(* * LANGUAGE : ANS Forth * PROJECT : Forth Environments * DESCRIPTION : Design of FIR filters. Mostly LPF. Uses remez(1) to optimize the response. * CATEGORY : EE Tools * AUTHOR : Marcel Hendrix * LAST CHANGE : Sunday, July 28, 2002 11:05 AM, Marcel Hendrix, Remez exchange algorithm * LAST CHANGE : July 14, 2002, Marcel Hendrix *) NEEDS -miscutil NEEDS -fsl_util NEEDS -sketcher REVISION -fir "ÄÄÄ Plot FIR output Version 1.10 ÄÄÄ" PRIVATES DOC (* From: 'Digital Filter Designer's Handbook', by C. Britton Rorabaugh This FIR filter tool knows about type 1..4 FIR's, here called Symmetric&Odd, Symmetric&Even, Asymmetric&Even, Asymmetric&Odd instead of 1, 2, 3 and 4. There are words to turn a supplied array into the discrete impulse response of an ideal lowpass, highpass, bandpass or bandstop filter. The number of filter taps (implied by the length of the supplied array) is unlimited. After you have the impulse response, another word allows to get the frequency response of that filter into some array. This Hd{ array can be plotted on the screen/printer/clipboard or dumped to a file with FirResponse and }FIRPlot. You can select linear (default) or logarithmic output with the word >Log. To improve the performance of the filter it is sometimes necessary to prefilter the data. This is done by embedding a convolution window into the filter's impulse response. The package has words that design an n-tap triangular, Hann, or Hamming window. The combined effect of an impulse response and a desired window can be studied with the FirWResponse word. DIRAC-PLOT and STEP-PLOT show the impulse and step response of a given FIR filter (wdirac-plot and wstep-plot allow to add a window first). You can try to match the above pre-cooked lowpass etc. designs using FreqSamplingFIR . This word allows to specify an arbitrary frequency characteristic and again builds a matching impulse response for you. The big problem with FreqSamplingFIR is that you have no control over the frequency response at points you _don't_ specify. This is analog to trying to fit an arbitrary polynomial through a limited number of points. A very nice word (well, it was very difficult to write :-) is therefore remez(1). For a type 1, i.e. Symmetric&Odd, FIR low-pass filter this word designs an equi-ripple Chebyshev-based filter with specifiable passband frequency, passband ripple, stopband frequency and stopband ripple. *) ENDDOC -- ------------------------------------------------------------------------------------------------------------ 0 VALUE /aptr PRIVATE 0 VALUE /bptr PRIVATE 0 VALUE /wptr PRIVATE DOUBLE DARRAY a1{ PRIVATE DOUBLE DARRAY b1{ PRIVATE DOUBLE DARRAY w1{ PRIVATE : OPEN-h(z) ( a{ b{ -- ) DUP CDIM TO /bptr SWAP DUP CDIM TO /aptr LOCALS| a{ b{ | /aptr /bptr MAX TO /wptr a1{ /wptr }malloc /aptr 0 ?DO a{ I } DF@ a1{ I } DF! LOOP b1{ /wptr }malloc /bptr 0 ?DO b{ I } DF@ b1{ I } DF! LOOP w1{ /wptr }malloc ; -- Evaluate the filter H(z) = b(z) / a(z) -- For a FIR filter a{ is zero. : h(z) ( F: xk -- yk ) 0e FLOCALS| yk xk | b1{ 0 } DF@ xk F* w1{ 0 } DF@ F+ FDUP TO yk /wptr 2- 0 ?DO b1{ I 1+ } DF@ xk F* a1{ I 1+ } DF@ yk F* F- w1{ I } DF+! LOOP b1{ /wptr 1- } DF@ xk F* a1{ /wptr 1- } DF@ yk F* F- w1{ /wptr 2- } DF! ; : CLOSE-h(z) ( -- ) a1{ }free b1{ }free w1{ }free CLEAR /aptr CLEAR /bptr CLEAR /wptr ; \ Simple IIR low-pass \ a{ 2 }RESIZE 1e a{ 0 } DF! -0.9e a{ 1 } DF! b{ 1 }RESIZE 1e b{ 0 } DF! -- ------------------------------------------------------------------------------------------------------------ 0 VALUE h1{ PRIVATE : OPEN-h(k) ( h{ -- ) DUP TO h1{ CDIM TO /wptr w1{ /wptr }malloc ; -- Evaluate hh : h(k) ( F: xk -- yk ) w1{ 0 } DF! 0e /wptr 0 ?DO w1{ /wptr 1- I - } DF@ h1{ I } DF@ F* F+ LOOP ( w0 ) /wptr 2 < ?EXIT 1 /wptr 1- DO w1{ I 1- } DF@ w1{ I } DF! -1 +LOOP ; : CLOSE-h(k) ( -- ) w1{ }free CLEAR h1{ CLEAR /wptr ; -- ------------------------------------------------------------------------------------------------------------ -- 0 <= lambda <= PI : >Lambda ( F: fcut fsample -- lambda ) F/ PI F* ; DEFER fir() ( F: lambda -- r ) ( hh{ -- ) : >Log ( mag{ -- ) DUP CDIM 0 ?DO DUP I } DUP DF@ FLOG 20e F* DF! LOOP DROP ; -- "Symmetric&Odd" means: h(nT) is symmetric around 0 and N is odd. :NONAME FDUP FLOCALS| lambda dlambda | DUP CDIM LOCALS| #taps hh{ | hh{ #taps 1- 2/ } DF@ 0 #taps 1- 2/ 1- DO hh{ I } DF@ F2* lambda FCOS F* F+ dlambda +TO lambda -1 +LOOP FABS ; =: Symmetric&Odd :NONAME FDUP F2/ FLOCALS| lambda dlambda | DUP CDIM LOCALS| #taps hh{ | 0e 0 #taps 2/ 1- DO hh{ I } DF@ F2* lambda FCOS F* F+ dlambda +TO lambda -1 +LOOP FABS ; =: Symmetric&Even :NONAME FDUP FLOCALS| lambda dlambda | DUP CDIM LOCALS| #taps hh{ | 0e 0 #taps 1- 2/ 1- DO hh{ I } DF@ F2* lambda FSIN F* F+ dlambda +TO lambda -1 +LOOP FABS ; =: AntiSymmetric&Odd :NONAME FDUP F2/ FLOCALS| lambda dlambda | DUP CDIM LOCALS| #taps hh{ | 0e 0 #taps 2/ 1- DO hh{ I } DF@ F2* lambda FSIN F* F+ dlambda +TO lambda -1 +LOOP FABS ; =: AntiSymmetric&Even Symmetric&Odd IS fir() -- Fill the output frequency array : FirResponse ( hh{ Hd{ -- ) DUP CDIM LOCALS| #points Hd{ hh{ | PI #points S>F F/ FLOCAL lambda #points 0 DO hh{ lambda I S>F F* fir() Hd{ I } DF! LOOP ; -- Define datawindows to suppress 'bleed' -------------------------------------------------------- -- A very bad window : MakeTriangularWindow ( window{ -- ) DUP CDIM DUP 1 AND LOCALS| odd? #taps window{ | odd? IF #taps 2/ 1+ 0 DO 1e I 2* S>F #taps 1+ S>F F/ F- FDUP window{ #taps 1- 2/ I + } DF! window{ #taps 1- 2/ I - } DF! LOOP ELSE #taps 2/ 0 DO 1e I 2* 1+ S>F #taps 1+ S>F F/ F- FDUP window{ #taps 2/ I + } DF! window{ #taps 2/ 1- I - } DF! LOOP ENDIF ; -- The Von Hann window is OK. : MakeHannWindow ( window{ -- ) DUP CDIM DUP 1 AND LOCALS| odd? #taps window{ | odd? IF #taps 2/ 1+ 0 DO I S>F PI*2 F* #taps 1- S>F F/ FCOS F2/ 0.5e F+ FDUP window{ #taps 1- 2/ I + } DF! window{ #taps 1- 2/ I - } DF! LOOP ELSE \ or 1- #taps 2/ 0 DO I 2* 1+ S>F PI F* #taps 1- S>F F/ FCOS F2/ 0.5e F+ FDUP window{ #taps 2/ I + } DF! window{ #taps 2/ 1- I - } DF! LOOP ENDIF ; -- The Hamming window is very good. : MakeHammingWindow ( window{ -- ) DUP CDIM DUP 1 AND LOCALS| odd? #taps window{ | odd? IF #taps 2/ 1+ 0 DO I S>F PI*2 F* #taps 1- S>F F/ FCOS 0.46e F* 0.54e F+ FDUP window{ #taps 1- 2/ I + } DF! window{ #taps 1- 2/ I - } DF! LOOP ELSE #taps 2/ 0 DO I 2* 1+ S>F PI F* #taps 1- S>F F/ FCOS 0.46e F* 0.54e F+ FDUP window{ #taps 2/ I + } DF! window{ #taps 2/ 1- I - } DF! LOOP ENDIF ; DOUBLE DARRAY hh*window{ PRIVATE -- Compute the effect of combined windowing + filtering : FirWResponse ( hh{ window{ Hd{ -- ) DUP CDIM 3 PICK CDIM LOCALS| #taps #points Hd{ window{ hh{ | PI #points S>F F/ FLOCAL lambda hh*window{ #taps }malloc #taps 0 DO hh{ I } DF@ window{ I } DF@ F* hh*window{ I } DF! LOOP #points 0 DO hh*window{ lambda I S>F F* fir() Hd{ I } DF! LOOP hh*window{ }free ; -- Compute Ideal Filters -------------------------------------------------------------------------- -- Do NOT use with AntiSymmetric&xx because A(0) = 0 : IdealLowpass ( hh{ -- ) ( F: lambdaU -- ) IS? fir() DUP AntiSymmetric&Odd = SWAP AntiSymmetric&Even = OR ABORT" Lowpass needs Symmetric&Odd or Symmetric&Even" DUP CDIM LOCALS| #taps hh{ | FLOCAL lambdaU #taps 0 DO I S>F #taps 1- S>F F2/ F- FDUP F0= IF FDROP lambdaU PI F/ ELSE FDUP lambdaU F* FSIN FSWAP PI F* F/ ENDIF hh{ I } DF! LOOP ; -- Do NOT use with Symmetric&Even or AntiSymmetric&Odd because A(pi) = 0 : IdealHighpass ( hh{ -- ) ( F: lambdaU -- ) IS? fir() DUP Symmetric&Odd <> SWAP AntiSymmetric&Even <> AND ABORT" Highpass needs Symmetric&Odd or AntiSymmetric&Even" DUP CDIM LOCALS| #taps hh{ | FNEGATE FLOCAL -lambdaL #taps 0 DO I S>F #taps 1- S>F F2/ F- FDUP F0= IF FDROP 1e -lambdaL PI F/ F+ ELSE FDUP -lambdaL F* FSIN FSWAP PI F* F/ ENDIF hh{ I } DF! LOOP ; : IdealBandpass ( hh{ -- ) ( F: lambdaL lambdaU -- ) DUP CDIM LOCALS| #taps hh{ | FSWAP FNEGATE FLOCALS| -lambdaL lambdaU | #taps 0 DO I S>F #taps 1- S>F F2/ F- FDUP F0= IF FDROP lambdaU -lambdaL F+ PI F/ ELSE FDUP lambdaU F* FSIN FOVER -lambdaL F* FSIN F+ FSWAP PI F* F/ ENDIF hh{ I } DF! LOOP ; -- Do NOT use with Symmetric&Even or AntiSymmetric&xx because A(0) and/or A(pi) = 0 : IdealBandstop ( hh{ -- ) ( F: lambdaL lambdaU -- ) IS? fir() Symmetric&Odd <> ABORT" Bandstop needs Symmetric&Odd" DUP CDIM LOCALS| #taps hh{ | FLOCALS| lambdaU lambdaL | #taps 0 DO I S>F #taps 1- S>F F2/ F- FDUP F0= IF FDROP lambdaL lambdaU F- PI F/ F1+ ELSE FDUP lambdaL F* FSIN FOVER lambdaU F* FSIN F- FSWAP PI F* F/ ENDIF hh{ I } DF! LOOP ; -- From amplitudes at a fixed grid 0 <= Fs <= 0.5, construct the corresponding -- FIR impulse response. Supports all 4 filter types : FreqSamplingFIR ( amplitudes{ impulse{ firtype -- ) 0e FLOCAL x >S DUP CDIM LOCALS| mm hn{ ee{ | CASE S> Symmetric&Odd OF mm 0 ?DO ee{ 0 } DF@ I mm 1- 2/ - S>F PI*2 F* mm S>F F/ TO x mm 1- 2/ 1+ 1 ?DO x I S>F F* FCOS ee{ I } DF@ F* F2* F+ LOOP mm S>F F/ hn{ I } DF! LOOP ENDOF Symmetric&Even OF mm 0 ?DO ee{ 0 } DF@ I mm 1- 2/ - S>F PI*2 F* mm S>F F/ TO x mm 2/ 1+ 1 ?DO x I S>F F* FCOS ee{ I } DF@ F* F2* F+ LOOP mm S>F F/ hn{ I } DF! LOOP ENDOF AntiSymmetric&Odd OF mm 0 ?DO 0e mm 1- 2/ I - S>F PI*2 F* mm S>F F/ TO x mm 2/ 1+ 1 ?DO x I S>F F* FSIN ee{ I } DF@ F* F2* F+ LOOP mm S>F F/ hn{ I } DF! LOOP ENDOF AntiSymmetric&Even OF mm 0 ?DO ee{ mm 2/ } DF@ mm 1- 2/ I - S>F PI F* FSIN F* I mm 1- 2/ - S>F PI*2 F* mm S>F F/ TO x mm 2/ 1+ 1 ?DO x I S>F F* FSIN ee{ I } DF@ F* F2* F+ LOOP mm S>F F/ hn{ I } DF! LOOP ENDOF ENDCASE ; -- REMEZ Exchange Algorithm ----------------------------------------------------------------------- -- Static variables 0.215e FVALUE FreqPass PRIVATE -- where passband ends 0.315e FVALUE FreqStop PRIVATE -- where stopband starts 25e FVALUE rr PRIVATE -- ripple ratio (deltapass / deltastop) = 0.025 / 0.001. 0 VALUE gridMax PRIVATE #16 VALUE gridDensity PRIVATE -- grid resolution when searching, might always OK. #13 VALUE #r PRIVATE -- (mm + 1) / 2 #25 VALUE mm PRIVATE -- Filter length ( odd ) 0 VALUE fextr{ PRIVATE -- output: extremal frequencies 0 VALUE hn{ PRIVATE -- output: impulse response 0 VALUE =S PRIVATE -- functions in passband 0 VALUE =P PRIVATE -- functions in stopband DOUBLE DARRAY ee{ PRIVATE -- error, difference between wanted and actual response DOUBLE DARRAY x{ PRIVATE -- extremal frequencies DOUBLE DARRAY beta{ PRIVATE DOUBLE DARRAY gamma{ PRIVATE INTEGER DARRAY iFF{ PRIVATE -- frequency grid -- the frequency response we want is the following: : desLpfResp ( F: f -- A ) 0e FSWAP FreqPass F<= IF 1e F+ ENDIF ;P -- weighting function (stopband ripple is more important) : WeightLp ( F: f -- A ) 1e FSWAP FreqPass F<= IF rr F/ ENDIF ;P 0e FVALUE gP PRIVATE 0e FVALUE incP PRIVATE 0e FVALUE incS PRIVATE 0e FVALUE absDelta PRIVATE -- tricky; see remezSearch. Keep at 0? : gridInit ( -- ) freqPass freqStop F- 0.5e F+ #r S>F F/ FreqPass FSWAP F/ FLOOR F>S TO =P =P gridDensity * S>F TO gP #r 1+ =P - TO =S freqPass gP F/ TO incP 0.5e freqStop F- =S 1- gridDensity * S>F F/ TO incS =P =S + 1- gridDensity * 1+ TO gridMax ;P : gridFreq ( gI -- ) ( F: -- f ) S>F FDUP gP F<= IF incP F* ELSE gP F1+ F- incS F* freqStop F+ ENDIF ;P : initRemezA ( -- ) -1e FLOCAL sign 0e 0e 0e FLOCALS| numer denom delta | #r 1+ 0 ?DO iFF{ I } @ gridFreq PI*2 F* FCOS x{ I } DF! LOOP #r 1+ 0 ?DO sign FNEGATE TO sign 1e #r 0 ?DO I J <> IF x{ J } DF@ x{ I } DF@ F- F/ ENDIF LOOP ( alpha ) FDUP beta{ I } DF! I #r <> IF x{ I } DF@ x{ #r } DF@ F- F/ ENDIF ( alpha ) FDUP iFF{ I } @ gridFreq FDUP desLpfResp FROT F* +TO numer ( alpha freq ) weightLp F/ sign F* +TO denom LOOP numer denom F/ TO delta -1e TO sign #r 0 ?DO sign FNEGATE TO sign iFF{ I } @ gridFreq ( freq ) delta FOVER weightLP F/ sign F* FSWAP desLpfResp FSWAP F- gamma{ I } DF! LOOP ( delta FABS TO absDelta ) ;P : computeRemezA ( F: contFreq -- num/denom ) PI*2 F* FCOS 0e 0e FLOCALS| denom numer xCont | #r 0 ?DO xCont x{ I } DF@ F- ( term) FDUP FABS 1e-12 F< IF FDROP gamma{ I } DF@ UNLOOP EXIT ENDIF beta{ I } DF@ FSWAP F/ FDUP +TO denom gamma{ I } DF@ F* +TO numer LOOP numer denom F/ ;P : remezError ( -- ) 0e FLOCAL freq initRemezA gridMax 0 ?DO I gridFreq TO freq freq WeightLP freq desLpfResp freq computeRemezA F- F* ee{ I } DF! LOOP ;P : remezSearch ( -- ) 0e 0e 0e 0e FLOCALS| =e-1 =e0 =e+1 smallestVal | 0 LOCAL kk \ seach for extremum at f=0 ee{ 0 } DF@ TO =e0 ee{ 1 } DF@ TO =e+1 =e0 F0> =e0 =e+1 F> AND =e0 FABS absDelta F>= AND =e0 F0< =e0 =e+1 F< AND =e0 FABS absDelta F>= AND OR IF 0 iFF{ kk } ! 1 +TO kk ENDIF \ seach for extrema in passband gP F>S 1 ?DO ee{ I 1+ } DF@ TO =e+1 ee{ I } DF@ TO =e0 ee{ I 1- } DF@ TO =e-1 =e0 =e-1 F>= =e0 =e+1 F> AND =e0 F0> AND =e0 =e-1 F<= =e0 =e+1 F< AND =e0 F0< AND OR IF I iFF{ kk } ! 1 +TO kk ENDIF LOOP \ pick up an extremal frequency at passband edge gP F>S iFF{ kk } ! 1 +TO kk \ pick up an extremal frequency at stopband edge gP F>S 1+ iFF{ kk } ! 1 +TO kk \ seach for extrema in stopband gridMax gP F>S 2+ ?DO ee{ I 1+ } DF@ TO =e+1 ee{ I } DF@ TO =e0 ee{ I 1- } DF@ TO =e-1 =e0 =e-1 F>= =e0 =e+1 F> AND =e0 F0> AND =e0 =e-1 F<= =e0 =e+1 F< AND =e0 F0< AND OR IF I iFF{ kk } ! 1 +TO kk ENDIF LOOP \ test for extremum af f=0.5 ee{ gridMax } DF@ TO =e0 ee{ gridMax 1- } DF@ TO =e-1 =e0 F0> =e0 =e-1 F> AND =e0 FABS absDelta F>= AND =e0 F0< =e0 =e-1 F< AND =e0 FABS absDelta F>= AND OR IF gridMax iFF{ kk } ! 1 +TO kk ENDIF \ find and remove superfluous extremal frequencies kk #r 1+ <= ?EXIT kk #r - 1 ?DO ee{ iFF{ 0 } @ } DF@ FABS TO smallestVal 0 >S ( index_smallest ) kk 1 ?DO ee{ iFF{ I } @ } DF@ FABS FDUP smallestVal F< IF TO smallestVal -S I >S ELSE FDROP ENDIF LOOP -1 +TO kk kk S> ?DO iFF{ I 1+ } @ iFF{ I } ! LOOP LOOP ;P : PassbandRipple ( F: -- max min ) ee{ iFF{ 0 } @ } DF@ FABS FDUP FLOCALS| biggestVal smallestVal | =P 1 ?DO ee{ iFF{ I } @ } DF@ FABS FDUP smallestVal FMIN TO smallestVal biggestVal FMAX TO biggestVal LOOP biggestVal smallestVal ;P : StopbandRipple ( F: -- max min ) ee{ iFF{ =P } @ } DF@ FABS FDUP FLOCALS| biggestVal smallestVal | =S 1 ?DO ee{ iFF{ I =P + } @ } DF@ FABS FDUP smallestVal FMIN TO smallestVal biggestVal FMAX TO biggestVal LOOP biggestVal smallestVal ;P : remezStop2 ( -- bool ) PassbandRipple FOVER FSWAP F- FSWAP F/ 0.01e F< StopbandRipple FOVER FSWAP F- FSWAP F/ 0.01e F< AND ;P : remezFinish ( -- ) #r 1+ 0 ?DO iFF{ I } @ gridFreq fextr{ I } DF! LOOP #r 0 ?DO I S>F mm S>F F/ computeRemezA ee{ I } DF! LOOP ( abuses of ee -- not an error value ) ee{ hn{ IS? fir() FreqSamplingFIR ;P -- deltapass / deltastop = rr', the obtained ripple ratio (specied is rr). -- This only works for type 1, but it can be MADE to work for the other 3 quite easily. -- See the book by Rorabaugh, or the article by Parks and McClellan (1973). -- How to use this to design a highpass/bandpass/stopband ? : remez(1) ( f{ h{ n -- ) ( F: freqPass freqStop rr -- deltapass deltastop ) TO mm TO hn{ TO fextr{ TO rr TO freqStop TO freqPass mm 1+ 2/ TO #r gridInit beta{ mm }malloc gamma{ mm }malloc x{ mm }malloc hn{ mm }malloc fextr{ #r 1+ }malloc iFF{ gridMax 1+ }malloc ee{ gridMax 1+ }malloc freqPass gP F2* F/ +TO freqPass ( Initial guess of extremal frequencies ) =P 0 ?DO I 1+ gridDensity * iFF{ I } ! LOOP =S 0 ?DO I gridDensity * 1+ gP F>S + iFF{ I =P + } ! LOOP ( Find optimal location of extremal frequencies ) #21 1 DO remezError remezSearch remezStop2 ?LEAVE LOOP PassbandRipple FDROP 0.0e WeightLP F/ StopbandRipple FDROP 0.5e WeightLP F/ remezFinish ( destroys ee{ which is needed for ripple computation ) ee{ }free x{ }free gamma{ }free beta{ }free iFF{ }free ; -- For LPF only; not sure if it would work for other than type 1. : APPROX-N ( F: freqPass freqStop deltaPass deltaStop -- ) ( -- n ) F* FLOG -10e F* 13e F- -FROT FSWAP F- 14.6e F* F/ F1+ FROUND F>S 1 OR ( make odd) ; -- ------------------------------------------------------------------------------------------------ 0 VALUE addr PRIVATE : [plot] ( ix -- ) ( F: -- r ) addr SWAP } DF@ ;P : }FIRPlot ( Hd{ -- ) TO addr 0 addr CDIM 1- ['] [plot] 'SKETCHA ; [DEFINED] testing [IF] #10 DOUBLE ARRAY a{ #10 DOUBLE ARRAY b{ #41 DOUBLE ARRAY hh{ #41 DOUBLE ARRAY Fextremal{ #41 DOUBLE ARRAY window{ #400 DOUBLE ARRAY Hd{ -- Test the quality of the four possible windows. Insightful! : WINDOWINGTEST ( -- ) PMULTI -90e 3e PSCALE White TO sketch-color hh{ Hd{ FirResponse Hd{ >Log Hd{ }FIRPlot Yellow TO sketch-color window{ MakeTriangularWindow hh{ window{ Hd{ FirWResponse Hd{ >Log Hd{ }FIRPlot Green TO sketch-color window{ MakeHannWindow hh{ window{ Hd{ FirWResponse Hd{ >Log Hd{ }FIRPlot Red TO sketch-color window{ MakeHammingWindow hh{ window{ Hd{ FirWResponse Hd{ >Log Hd{ }FIRPlot PSINGLE ;P -- Example 1: Always SYMMETRIC; n-tap FIR lowpass : LPF-TEST ( -- ) Symmetric&Odd IS fir() hh{ #21 }RESIZE window{ #21 }RESIZE 2000e 5000e >Lambda hh{ IdealLowpass WINDOWINGTEST ; -- Example 2: SYMMETRIC and ODD or ANTISYMMETRIC and EVEN; n-tap FIR highpass : HPF-TEST ( -- ) Symmetric&Odd IS fir() hh{ #21 }RESIZE window{ #21 }RESIZE 3000e 5000e >Lambda hh{ IdealHighpass WINDOWINGTEST ; -- Example 3: No restrictions; n-tap FIR bandpass : BPF-TEST ( -- ) Symmetric&Odd IS fir() hh{ #21 }RESIZE window{ #21 }RESIZE 2000e 5000e >Lambda 3000e 5000e >Lambda hh{ IdealBandpass WINDOWINGTEST ; -- Example 4: always ODD; n-tap FIR bandstop : BSF-TEST ( -- ) Symmetric&Odd IS fir() hh{ #31 }RESIZE window{ #31 }RESIZE 2000e 5000e >Lambda 3000e 5000e >Lambda hh{ IdealBandstop WINDOWINGTEST ; #400 DOUBLE ARRAY impulse{ PRIVATE : DIRAC-PLOT ( -- ) hh{ OPEN-h(k) 1e h(k) impulse{ 0 } DF! hh{ CDIM 1 ?DO 0e h(k) impulse{ I } DF! LOOP impulse{ TO addr 0 hh{ CDIM 1- ['] [plot] 'SKETCHA CLOSE-h(k) ; : STEP-PLOT ( -- ) hh{ OPEN-h(k) hh{ CDIM 0 ?DO 1e h(k) impulse{ I } DF! LOOP impulse{ TO addr 0 hh{ CDIM 1- ['] [plot] 'SKETCHA CLOSE-h(k) ; : WDIRAC-PLOT ( -- ) hh*window{ hh{ CDIM }malloc hh{ CDIM 0 ?DO hh{ I } DF@ window{ I } DF@ F* hh*window{ I } DF! LOOP hh*window{ OPEN-h(k) 1e h(k) impulse{ 0 } DF! hh*window{ CDIM 1 ?DO 0e h(k) impulse{ I } DF! LOOP impulse{ TO addr 0 hh*window{ CDIM 1- ['] [plot] 'SKETCHA hh*window{ }free CLOSE-h(k) ; : WSTEP-PLOT ( -- ) hh*window{ hh{ CDIM }malloc hh{ CDIM 0 ?DO hh{ I } DF@ window{ I } DF@ F* hh*window{ I } DF! LOOP hh*window{ OPEN-h(k) hh*window{ CDIM 0 ?DO 1e h(k) impulse{ I } DF! LOOP impulse{ TO addr 0 hh*window{ CDIM 1- ['] [plot] 'SKETCHA hh*window{ }free CLOSE-h(k) ; [THEN] :ABOUT CR ." Usage: Symmetric&Odd | Symmetric&Even | Asymmetric&Even | Asymmetric&Odd IS fir()" CR ." ( window{ -- ) MakeTriangularWindow | MakeHannWindow | MakeHammingWindow" CR ." ( hh{ f_lambda -- ) IdealLowpass | IdealHighpass -- design filter in hh{" CR ." ( hh{ f_lambdaL f_lambdaU -- ) IdealBandpass | IdealBandstop -- design filter in hh{" CR ." ( F: fcut fsample -- lambda ) >Lambda -- normalized filter corner frequencies" CR ." ( hh{ Hd{ -- ) FirResponse -- gain without windowing; fir() must be set" CR ." ( hh{ window{ Hd{ -- ) FirWResponse -- gain with windowing; fir() must be set" CR ." ( Hd{ h{ firtype -- ) FreqSamplingFIR -- design FIR impulse response from A(f)" CR ." ( h{ -- ) OPEN-h(k) -- prepare for h(k)" CR ." ( F: xk -- yk ) h(k) -- evaluate h(k) in time domain" CR ." CLOSE-h(k) -- stop using h(k)" CR ." ( a{ b{ -- ) OPEN-h(z) -- prepare for h(z)" CR ." ( F: xk -- yk ) h(z) -- evaluate h(z) in time domain" CR ." CLOSE-h(z) -- stop using h(z)" CR ." ( Hd{ -- ) >Log -- modifies Hd{ permanently" CR ." ( Hd{ -- ) }FIRPlot -- plot amplitude response" CR ." ( f{ h{ n -- ) -- odd-length (n) LPF FIR, ripple ratio rr" CR ." ( F: fP fS rr -- dP dS ) remez(1) -- size f{ = n+1/2, h{ = n, f{ is extreme freqs." CR ." -- fir() should be correctly set up for type 1" CR ." Note: hh{ is the sampled impulse response, Hd{ the computed frequency response" CR CR ." Available when 1 =: testing is set before loading..." CR ." ( arrays ) a{ b{ hh{ Hd{ window{ Fextremal{" CR ." Demos: LPF-TEST | HPF-TEST | BPF-TEST | BSF-TEST -- check effect of windowing on these filters" CR ." dirac-plot | step-plot -- after hh{ is valid" CR ." wdirac-plot | wstep-plot -- after hh{ and window{ are valid" ; DEPRIVE .ABOUT -fir CR (* End of Source *)