On Thu, 9 Nov 2000, Lawrence Lile wrote:

> > I guess it depends on your application... Which is?
> >
> >
> > Scott
>
>
>
> I'm attempting to pick out ~~20KHZ audio tones from the  random sounds that
> are bouncing around a room just above our hearing limits.  I've got a
> control that detects [patentable subject matter deleted] based on the
> presence of absence of an audio tone.  I want to do two things to make sure
> that the tone I detect isn't just noise:  1st bandpass filter it so I know I
> have the correct frequency, and then  also gate it so I know it arrived at
> the right time, just after I sent the tone out of a speaker.

20kHz audio? You've got good ears! (My hearing falls off 60dB at 2kHz in one ear
and 60dB at 3kHz in the other.)

>
> I need to do this on the cheap which is why I am even considering a digital
> filter.  I haven't tried yet to figure out if I can process information fast
> enough in a lowly PIC (and I'll be using the lowliest of all if I can..)
>
> I am still trying to figure out something.  Was your [ever so useful]
> diatribe on DFT recommending detecting a sine wave by comparing with square
> waves of the desired frequency?  or detecting square waves by comparing with
> square waves?  Sorta like detecting round pegs with square holes.....

The algorithm attempts to detect sine waves - although it works for square waves
as well (in fact probably even better). The idea is that a sine wave can be
converted into a square wave by passing it through a comparator. If the
comparator switches states at the zero crossing, then a sine wave without any
noise would be converted into a 50% duty cycle square wave. Now if you recall
the Fourier Series for a square wave (see
http://www.dattalo.com/technical/theory/sqwave.html ) with 50% duty cycle is
absent of even harmonics and the odd harmonic strength is inversely proportional
to the harmonic number (e.g. the third harmonic is 1/3 as strong as the
fundamental, etc.). As you can probably guess, these harmonics complicate the
detection algorithm. The is one of the reasons why in the one application I
mentioned I separated the processing of the two tones.

 >
> It seems like what you are implying in the following code snippet is that
> the input data has simply been run through a comparator so the input sine
> wave is now just a square wave of the same frequency..... I'll have to model
> this in my spreadsheet and see how it works out.
>
> hmmm... a comparator is one of the  the CHEAPEST A/D converters ever
> conceived....   1/6 hex inverter at $0.15/6 each!
>
>
> Scott wrote:
> >In PIC parlance,
>
>
> >sum_of_products    equ  0x20   ;A register variable
> >input_sample       equ  0x21   ;In LS bit
> >square_wave_sample equ  0x22   ;In LS bit
>
> >        movf    input_sample,W            <----------   This is just a 1 or
> 0, no?
> >        xorwf   square_wave_sample,W
> >        skpnz
> >         movlw  -2
> >        addlw   1
> >        addwf   sum_of_products,F

Well, that is the slow way to do it :). If you're trying to detect just one tone
then there's room for optimization. Recall that in addition to computing the DFT
by counting the pulses you also need to synthesize the quadrature wave forms. If
you observe the quadrature waveforms:

  Q3 Q0 Q1 Q2 Q3 Q0 Q1 Q2 Q3 Q0 Q1 Q2 Q3 Q0 Q1 Q2 Q3
    +-----+     +-----+     +-----+     +-----+     +-----+
C   |     |     |     |     |     |     |     |     |     |
  --+     +-----+     +-----+     +-----+     +-----+

       +-----+     +-----+     +-----+     +-----+     +-----+
S      |     |     |     |     |     |     |     |     |     |
  -----+     +-----+     +-----+     +-----+     +-----+

You'll notice that there are four cases: Both high, Both low, first high second
low and the second high with the first low. Now, If you created an isochronous
loop that takes exactly one period of your tone to execute, then you could break
the loop into these four cases. I labeled these 4 states, or quadratures as Q0
through Q3:

Q0  - C is high S is low
Q1  - both high
Q2  - C is low S is high
Q3  - both low

In your specific case, f=20kHz, the time for one cycle is 50uS. Suppose you
choose the PIC oscillator to be a 16MHz crystal. The instruction cycle time is
250nS. The number of instructions that can execute in 50uS is 200 or 12.5uS = 50
instructions for each quadrature.

Recall also that the `dot product' is simply:

      N
    ----
    \   f * q
    /
    ----

Where f is the unknown signal were trying to measure and q is one of the (known)
quadrature wave forms. Now f is analog signal that has been quantized by the
comparator. So what the pic will see is a train of pulses or a square wave that
has the same frequency of the 20kHz sine wave. As I pointed out on my web page,
the multiplication boils down to a simple boolean AND.

In psuedo code, you may write it like so:


  C_dot_product = 0;
  S_dot_product = 0;

  for(i=0; i<NUMBER_OF_SAMPLES_PER_CYCLE; i++) {

    if(C is high AND f is high)
      C_dot_product++;

    if(S is high AND f is high)
      S_dot_product++;

    synthesize next sample for C and S
  }

Or if we apply the optimization of breaking the dot product into 4 pieces:


do {

  Q0_dp = 0;
  Q1_dp = 0;
  Q2_dp = 0;

  for(i=0; i<NUMBER_OF_SAMPLES_PER_CYCLE/4; i++) {
    if(f is high)
      Q0_dp++;
  }

  for(i=0; i<NUMBER_OF_SAMPLES_PER_CYCLE/4; i++) {
    if(f is high)
      Q1_dp++;
  }

  for(i=0; i<NUMBER_OF_SAMPLES_PER_CYCLE/4; i++) {
    if(f is high)
      Q2_dp++;
  }

  for(i=0; i<NUMBER_OF_SAMPLES_PER_CYCLE/4; i++) {
    // don't need to do anything here, because the
    // dot product is always zero.
    // So, the quad is where we normalize the results
    // from the other three quads.
    C_dot_product += (Q0_dp + Q1_dp);
    S_dot_product += (Q1_dp + Q2_dp);

    Normalize();   // as described on the web page.

  }

} while(1);

Notice that the C and S square waves no longer need to be synthesized. Their
states are implied by the fact of where/when the code is excuting.

Of course, you wouldn't go off and write this in C, right? So in assembly, the
fastest thing to do would be to unwind the loops and write this:

;written off the top of my head (UNTESTED!!!)

    clrf   Q0_dp    ; Clear the initial dot products
    clrf   Q1_dp
    clrf   Q2_dp
    clrf   C_dp_lo
    clrf   C_dp_hi
    clrf   S_dp_lo
    clrf   S_dp_hi


; The macro for forming a dot_product over one quadrature
; Note that there are 2-instructions per sample. And since
; the pic-crystal is 16Mhz, we have 50us per quadrature or
; 200 instructions per quadrature. This gives 100 samples

#define SAMPLES_PER_QUAD  100

MACRO  Q_DP  q

     variable i
i = 0

     while i < SAMPLES_PER_QUAD
      BTFSC  IOPORT,IOBIT
       INCF  q,f
i += 1
     endw

  endm



Q0:

    Q_DP Q0_dp

Q1:
    Q_DP Q2_dp

Q2:
    Q_DP Q1_dp

Q3:

;  C_dp += (Q0_dp + Q1_dp)

    movf  Q0_dp,w
    addwf Q1_dp,w
    skpnc
     incf C_dp_hi,f
    addwf C_dp_lo,f
    skpnc
     incf C_dp_hi,f

;  S_dp += (Q1_dp + Q2_dp)

    movf  Q1_dp,w
    addwf Q2_dp,w
    skpnc
     incf S_dp_hi,f
    addwf S_dp_lo,f
    skpnc
     incf S_dp_hi,f

 ...

  normalize...

  check if we should quit...

    goto  Q0


I'll save the normalization for another post...

The Q_DP macro is invoked for the first three quadratures. This will sample the
input 100 times and increment the Qx_dp counters. For the last quadrature, we
take advantage of the fact that both of the synthesized square waves are low and
thus proceed to normalize the results. Part of the normalization should also
include a recursive low pass filter like those recently discussed to smooth out
the algorithm over many 20kHz cycles. A 20kHz waveform is declared detected when
the filtered normalized result exceeds an emperically determined threshold.

Is this clear?

Scott

--
http://www.piclist.com hint: The PICList is archived three different
ways.  See http://www.piclist.com/#archives for details.