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 )