            Path: typhoon.bart.nl!newshub.vianetworks.nl!newshub.vianetworks.de!newshub.bart.net!213.51.129.3.MISMATCH!newshub1.home.nl!home.nl!news.cambrium.nl!news.cambrium.nl!news2.euro.net!newsfeed.online.be!sn-xit-02!sn-xit-01!sn-post-01!supernews.com!news.supernews.com!not-for-mail
            From: andrew29@littlepinkcloud.invalid
            Newsgroups: comp.lang.forth
            Subject: Re: Application Development
            Date: Wed, 31 Mar 2004 20:13:47 -0000
            Message-ID: <106m9jrfi4ccu65@news.supernews.com>
            Sender: Andrew Haley <aph@cuddles.cambridge.redhat.com>
            References: <4061561b_3@127.0.0.1>   <c440jq$2e30ko$1@ID-217727.news.uni-berlin.de> <c441dr$3uu$07$1@news.t-online.com> <4065d8a7$0$3057$61fed72c@news.rcn.com> <6vejj1-gpf.ln1@cohen.paysan.nom> <1gbezjl.7n9acbqzxayuN%ward@megawolf.com> <4068e303$0$3069$61fed72c@news.rcn.com> <1gbglv0.1d0rmra14xtfk8N%ward@megawolf.com> <106jg8j4bepe508@news.supernews.com> <1gbgrzg.1b6kpce9mkvi8N%ward@megawolf.com> <106jjbpac3h7afd@news.supernews.com> <1gbgvg3.3eb1zsr3jd50N%ward@megawolf.com>
            User-Agent: tin/1.4.4-20000803 ("Vet for the Insane") (UNIX) (Linux/2.4.20-20.9 (i686))
            X-Complaints-To: abuse@supernews.com
            Lines: 154
            X-Received-Date: Wed, 31 Mar 2004 22:15:14 MET DST (typhoon.bart.nl)
            Xref: typhoon.bart.nl comp.lang.forth:46777

            [..] I've attached a little Forth program that shows how to do
            resampling correctly.  Its output with a sample stream at 0.95 *
            the Nyquist frequency looks like this:

                          #
                                 #
                                        #
                                               #
                                                    *  .809016994383555
                                                       #
                                                       #
                                                    #
                                                #
                                         #
                                  #
                           #
                     *  -.707106781197038
                  #
                #
                 #
                     #
                           #
                                  #
                                         #
                                               *  .587785252304614
                                                    #
                                                      #
                                                      #
                                                    #
                                               #
                                        #
                                 #
                          *  -.453990499753049
                     #
                  #
                 #
                  #
                      #
                            #
                                   #
                                          *  .309016994389529
                                                #
                                                    #
                                                      #
                                                      #
                                                   #
                                              #
                                       #
                                *  -.156434465055553
                          #
                     #
                 #
                 #
                   #
                       #
                             #

            The sample points are marked with their numeric values, and the
            intermediate values are interpolated by the filter, which is a 128
            point minimum phase FIR filter with 8 times oversampling.  This is
            rather crude in comparison with the real world versions used in
            commercial DACs, but it works well enough.  It gets rather ragged as
            you approach the Nyquist frequency, but no filter is perfect.

            The sampled data is in an array called samples.  The word init
            generates a sampled sine wave at 0.95 the Nyquist frequency, but you
            can put anything in the samples array and see what the reconstruction
            filter does.  If you put random data in there you'll get out white
            noise low pass filtered at the Nyquist frequency.

            Andrew.

    : FARRAY ( n -- )
            CREATE  FLOATS ALLOT IMMEDIATE
            DOES>   POSTPONE LITERAL EVAL" []FLOAT " ;

    :INLINE sinc ( F: f -- f )
            FDUP  F0= IF  FDROP 1e
                    ELSE  PI F*  FDUP FSIN  FSWAP F/
                   ENDIF ;

    #64 =: windowlen
      8 =: oversampling

    windowlen 2* #32 + =: datasize
    datasize FARRAY samples

    :INLINE window ( F: f -- f )
            FABS  FDUP F>S
            windowlen > IF  0e
                      ELSE  windowlen S>F  F/  1e FSWAP F-
                     ENDIF ;

    :INLINE tap ( F: f -- f )  FDUP sinc FSWAP window F* ;

    : init  ( -- )
            datasize 0 DO   I S>F 0.95e F*  PI F*  FSIN
                            I samples F!
                     LOOP ;

            init

    :INLINE ?int ( F: f -- t )  FDUP F>S S>F F= ;

    : show ( F: ft fy -- )
            FDUP  1.2e F+  20e F*
            F>S SPACES
            FSWAP ?INT IF  ." *  " F.
                     ELSE  ." # "  FDROP
                    ENDIF
            CR ;

    : filter ( F: ft -- f )
            0e FLOCAL sum
            FDUP F>S  DUP windowlen +  SWAP windowlen -
                DO  FDUP I S>F F- tap
                    I samples F@ F*  +TO sum
              LOOP
            FDROP  sum ;

    : demo ( -- )
            datasize windowlen - oversampling *
            windowlen oversampling *
              DO  I S>F  oversampling S>F F/
                  FDUP filter show
            LOOP ;

            CR .( Try: show )