Loading ...
Sorry, an error occurred while loading the content.

Re: [buildcheapeeg] need help with DSP

Expand Messages
  • Jim Peters
    ... Here we go ... When you do an FFT of a block of sound, you can t just FFT it straight off. Well, you can, but you end up with lots of distortion because
    Message 1 of 11 , Jun 15, 2002
    • 0 Attachment
      Michal Wallace wrote:
      > What I tried to do was create 10 seconds of fake data,
      > simulating "pure" brainwaves (such as a steady 6Hz for
      > theta), as well as some contrived mixtures. Then, I
      > ran an FFT on them to get the spectrogram.
      >
      > The output bears some resemblance to what I wanted, but
      > it's grossly distorted, and I'm not sure what's wrong.
      >
      > I think the code *MIGHT* be working fine except that it's
      > using the wrong scale for the "bins", because it looks like
      > it might be right, except it's all squished into the top
      > bars and there seems to be a lot of noise.

      Here we go ...

      When you do an FFT of a block of sound, you can't just FFT it straight
      off. Well, you can, but you end up with lots of distortion because
      you have effectively cut a rectangle out of the sound, and there are
      big clicks at the start and end of it. (And when you give a chunk of
      sound to an FFT, it acts as if that sound chunk is looped forever, so
      you can imagine what kind of distortion that creates if you don't tidy
      up the edges).

      So, you need to apply a 'window' to the data before FFTing it. If you
      look at BWView, if you click somewhere on the display, you can see the
      'window function' being used in the display above. This looks
      something like a bell-shape.

      So, if you want nice results from your FFT, you need to window the
      data first -- which means multiplying it by a 'window function',
      effectively flattening off the edges. Don't think you can choose just
      any old function, because each different shape has its own particular
      frequency response. There are several types of window, but I used the
      'Blackman' window in my code.

      You mention that you've already looked at the code, so you should have
      no trouble finding this (in analysis.c around line 410 onwards). The
      multiplier is given by this:

      double mag= 0.42 + 0.5 * cos(ang) + 0.08 * cos(2*ang); // Blackman window

      Where 'ang' goes from -pi at one end of the chunk to +pi at the other.
      It's a while since I've coded in Python, so I'll prototype it in C for
      you:

      for (a= 0; a<64; a++) {
      double ang= (a + 0.5 - 32.0) * M_PI / 32.0;
      double mag= 0.42 + 0.5 * cos(ang) + 0.08 * cos(2*ang);
      slice[a] *= mag;
      }

      I think that should be correct.

      Note that a 6Hz wave is not going to be very easy to see if you are
      only using a 0.25 second window (that's only 1.5 of a wavelength
      you're expecting it to recognise!). Bear in mind that the lowest band
      in your FFT will be centred on 0Hz, then 4Hz, then 8Hz, 16Hz, 24Hz,
      and so on up to 128Hz.

      Jim

      --
      Jim Peters (_)/=\~/_(_) jim@...
      (_) /=\ ~/_ (_)
      Uazú (_) /=\ ~/_ (_) http://
      B'ham, UK (_) ____ /=\ ____ ~/_ ____ (_) uazu.net
    • Jim Peters
      ... Oops -- my mistake. It goes 0, 4, 8, 12, 16, 20, ... 128; basically (0..32)*4Hz. The 4Hz comes from the window-length: 256/64 == 4Hz. I hope I haven t
      Message 2 of 11 , Jun 15, 2002
      • 0 Attachment
        Jim Peters wrote:
        > Bear in mind that the lowest band in your FFT will be centred on
        > 0Hz, then 4Hz, then 8Hz, 16Hz, 24Hz, and so on up to 128Hz.

        Oops -- my mistake. It goes 0, 4, 8, 12, 16, 20, ... 128; basically
        (0..32)*4Hz. The 4Hz comes from the window-length: 256/64 == 4Hz.

        I hope I haven't made any other slips --

        Jim

        --
        Jim Peters (_)/=\~/_(_) jim@...
        (_) /=\ ~/_ (_)
        Uazú (_) /=\ ~/_ (_) http://
        B'ham, UK (_) ____ /=\ ____ ~/_ ____ (_) uazu.net
      • Michal Wallace
        ... Aha! I wondered why one of the samples came out clean and the others didn t.. It s because that one, when looped, was a continuous wave and the others had
        Message 3 of 11 , Jun 15, 2002
        • 0 Attachment
          On Sat, 15 Jun 2002, Jim Peters wrote:

          > When you do an FFT of a block of sound, you can't just FFT it straight
          > off. Well, you can, but you end up with lots of distortion because
          > you have effectively cut a rectangle out of the sound, and there are
          > big clicks at the start and end of it. (And when you give a chunk of
          > sound to an FFT, it acts as if that sound chunk is looped forever, so
          > you can imagine what kind of distortion that creates if you don't tidy
          > up the edges).


          Aha! I wondered why one of the samples came out clean and
          the others didn't.. It's because that one, when looped, was
          a continuous wave and the others had the edges cut
          off. Thanks! :)

          > So, if you want nice results from your FFT, you need to
          > window the data first -- which means multiplying it by a
          > 'window function', effectively flattening off the edges.
          > Don't think you can choose just any old function, because
          > each different shape has its own particular frequency
          > response. There are several types of window, but I used
          > the 'Blackman' window in my code.

          Thanks. I'd seen the term window used in the various
          articles I read, but I hadn't understood what it was for
          until you just explained it.

          I took a look at your code, and then through the docs for
          Numpy, and found that it has an "MLab" (Matlab-compatible?)
          module which generates it:

          import MLab
          window = MLab.blackman(length)

          Then I can do "slice * window" in the loop and not have
          to recalculate the window each time.

          I tried it out, and it definitely cleans up the noise.
          Thanks! :)


          Cheers,

          - Michal http://www.sabren.net/ sabren@...
          ------------------------------------------------------------
          Switch to Cornerhost! http://www.cornerhost.com/
          High Powered Hosting - With a Human Touch. :)
          ------------------------------------------------------------
        • Michal Wallace
          ... Hey Jim, Thanks for this, too. I think I m starting to understand the properties of the FFT results a little better. Let me see if I have this straight: 1.
          Message 4 of 11 , Jun 15, 2002
          • 0 Attachment
            On Sat, 15 Jun 2002, Jim Peters wrote:

            > Jim Peters wrote:
            > > Bear in mind that the lowest band in your FFT will be centred on
            > > 0Hz, then 4Hz, then 8Hz, 16Hz, 24Hz, and so on up to 128Hz.
            >
            > Oops -- my mistake. It goes 0, 4, 8, 12, 16, 20, ... 128;
            > basically (0..32)*4Hz. The 4Hz comes from the
            > window-length: 256/64 == 4Hz.


            Hey Jim,

            Thanks for this, too. I think I'm starting to understand the
            properties of the FFT results a little better.

            Let me see if I have this straight:

            1. We want to pick up waves in the 0-Beta range, where
            beta is around 30. Some people may even want "super beta",
            so we may want to go up as high as 40-50?

            2. Since Nyquist's rule says we need to sample at least
            twice as fast as the highest frequency we want to detect,
            the hardware folks have chosen a sample rate of 128Hz
            or 256Hz. (I'm not sure about this, I thought I'd read
            it somewhere... Can anyone tell me the actual rates?)

            3. The FFT returns an array of the same length as the
            input, but we discard half of it because it deals with
            imaginary numbers. Each slot in the array (or "bin", I've
            seen it called) represents a set of frequencies
            max-detectible-frequency/window-length wide.
            Is that true?

            So, if I want each bar to represent 0hz, 1hz, ... 32hz, and
            we have 256 samples per second, then then I need to have all
            256 samples, which gives me 128 bins, each 1Hz wide. Then I
            can just show the first 32. Is that right?

            If so... Do I just update the screen once a second? Or
            should I show the last second's worth of data every 1/4th of
            a second, even though 75% of the data is the same each time?


            Thanks,

            - Michal http://www.sabren.net/ sabren@...
            ------------------------------------------------------------
            Switch to Cornerhost! http://www.cornerhost.com/
            High Powered Hosting - With a Human Touch. :)
            ------------------------------------------------------------
          • Jim Peters
            ... It would be better to make it flexible -- it might be 128Hz, 200Hz, or 256Hz -- who knows ? ... Yes, this all sounds correct to me so far. Remember that
            Message 5 of 11 , Jun 15, 2002
            • 0 Attachment
              Michal Wallace wrote:
              > 2. Since Nyquist's rule says we need to sample at least
              > twice as fast as the highest frequency we want to detect,
              > the hardware folks have chosen a sample rate of 128Hz
              > or 256Hz. (I'm not sure about this, I thought I'd read
              > it somewhere... Can anyone tell me the actual rates?)

              It would be better to make it flexible -- it might be 128Hz, 200Hz, or
              256Hz -- who knows ?

              > 3. The FFT returns an array of the same length as the
              > input, but we discard half of it because it deals with
              > imaginary numbers. Each slot in the array (or "bin", I've
              > seen it called) represents a set of frequencies
              > max-detectible-frequency/window-length wide.
              > Is that true?
              >
              > So, if I want each bar to represent 0hz, 1hz, ... 32hz, and
              > we have 256 samples per second, then then I need to have all
              > 256 samples, which gives me 128 bins, each 1Hz wide. Then I
              > can just show the first 32. Is that right?

              Yes, this all sounds correct to me so far.

              Remember that you are always trading frequency resolution against
              speed of response. You could have very narrow bins, but you would
              need a huge long window; or alternatively, you could have a very fast
              response, but you would end up with big fat bins.

              > If so... Do I just update the screen once a second? Or
              > should I show the last second's worth of data every 1/4th of
              > a second, even though 75% of the data is the same each time?

              It sounds like quite a good idea to show a second's worth of data 4
              times a second (or more), as you suggest.

              Jim

              --
              Jim Peters (_)/=\~/_(_) jim@...
              (_) /=\ ~/_ (_)
              Uazú (_) /=\ ~/_ (_) http://
              B'ham, UK (_) ____ /=\ ____ ~/_ ____ (_) uazu.net
            • larryjanow
              Hi Michal, You can t discard the imaginary half of the FFT output, you need it to calculate the actual output amplitude at each frequency. The formula is:
              Message 6 of 11 , Jun 17, 2002
              • 0 Attachment
                Hi Michal,

                You can't discard the 'imaginary' half of the FFT output, you need
                it to calculate the actual output amplitude at each frequency.

                The formula is:

                Amplitude(m) = square root of( ( Real(m) * Real(m) ) + ( Imag(m) *
                Imag(m) ) ) / ( N / 2 ), where N is the FFT sample size (256 in your
                example below)

                For example, to calculate the 10Hz amplitude using your example
                below, Real(m) would be bin 11 (bin 0 is the real part of the DC
                component), and Imag(m) would be bin 139 (bin 128 is the imaginary
                part of DC).

                You can also use the imaginary bins to calculate the phase of each
                component:

                phase = inverse tangent of( Imag(m) / Real(m) )

                Good luck!

                - Larry


                --- In buildcheapeeg@y..., Michal Wallace <sabren@m...> wrote:
                > On Sat, 15 Jun 2002, Jim Peters wrote:
                >
                > > Jim Peters wrote:
                > > > Bear in mind that the lowest band in your FFT will be centred
                on
                > > > 0Hz, then 4Hz, then 8Hz, 16Hz, 24Hz, and so on up to 128Hz.
                > >
                > > Oops -- my mistake. It goes 0, 4, 8, 12, 16, 20, ... 128;
                > > basically (0..32)*4Hz. The 4Hz comes from the
                > > window-length: 256/64 == 4Hz.
                >
                >
                > Hey Jim,
                >
                > Thanks for this, too. I think I'm starting to understand the
                > properties of the FFT results a little better.
                >
                > Let me see if I have this straight:
                >
                > 1. We want to pick up waves in the 0-Beta range, where
                > beta is around 30. Some people may even want "super beta",
                > so we may want to go up as high as 40-50?
                >
                > 2. Since Nyquist's rule says we need to sample at least
                > twice as fast as the highest frequency we want to detect,
                > the hardware folks have chosen a sample rate of 128Hz
                > or 256Hz. (I'm not sure about this, I thought I'd read
                > it somewhere... Can anyone tell me the actual rates?)
                >
                > 3. The FFT returns an array of the same length as the
                > input, but we discard half of it because it deals with
                > imaginary numbers. Each slot in the array (or "bin", I've
                > seen it called) represents a set of frequencies
                > max-detectible-frequency/window-length wide.
                > Is that true?
                >
                > So, if I want each bar to represent 0hz, 1hz, ... 32hz, and
                > we have 256 samples per second, then then I need to have all
                > 256 samples, which gives me 128 bins, each 1Hz wide. Then I
                > can just show the first 32. Is that right?
                >
                > If so... Do I just update the screen once a second? Or
                > should I show the last second's worth of data every 1/4th of
                > a second, even though 75% of the data is the same each time?
                >
                >
                > Thanks,
                >
                > - Michal http://www.sabren.net/ sabren@m...
                > ------------------------------------------------------------
                > Switch to Cornerhost! http://www.cornerhost.com/
                > High Powered Hosting - With a Human Touch. :)
                > ------------------------------------------------------------
              • larryjanow
                OOPS... The example I gave below would calculate the 11Hz amplitude, not 10Hz.
                Message 7 of 11 , Jun 17, 2002
                • 0 Attachment
                  OOPS...

                  The example I gave below would calculate the 11Hz amplitude, not
                  10Hz.

                  --- In buildcheapeeg@y..., "larryjanow" <larrywj@a...> wrote:
                  >
                  > For example, to calculate the 10Hz amplitude using your example
                  > below, Real(m) would be bin 11 (bin 0 is the real part of the DC
                  > component), and Imag(m) would be bin 139 (bin 128 is the imaginary
                  > part of DC).
                  >
                • Jim Peters
                  ... His python code looks like this: def makeSpectrogram(slice): assert len(slice)==64, we want 32 bins, so we need 64 samples res =
                  Message 8 of 11 , Jun 17, 2002
                  • 0 Attachment
                    larryjanow wrote:
                    > You can't discard the 'imaginary' half of the FFT output, you need
                    > it to calculate the actual output amplitude at each frequency.

                    His python code looks like this:

                    def makeSpectrogram(slice):
                    assert len(slice)==64, "we want 32 bins, so we need 64 samples"

                    res = abs(FFT.real_fft(slice))[:-1] # discard 33rd slot (is this okay?)
                    res = Numeric.floor(res) # round off to integers

                    assert len(res)==32, len(res)
                    return res

                    If real_fft is returning 33 values from a 64-sample input, then it
                    seems to me that it must be returning complex numbers. If so, then
                    the abs() function should automatically do the square root you
                    mention (i.e. sqrt(r*r+i*i)).

                    Perhaps this needs checking. However, obviously the Python
                    FFT.real_fft is not returning the imaginary parts as real numbers in
                    the places you are expecting it to, so your method won't work in this
                    particular case, and is based on some other FFT algorithm's output
                    layout.

                    Jim

                    --
                    Jim Peters (_)/=\~/_(_) jim@...
                    (_) /=\ ~/_ (_)
                    Uazú (_) /=\ ~/_ (_) http://
                    B'ham, UK (_) ____ /=\ ____ ~/_ ____ (_) uazu.net
                  • larryjanow
                    I made my comments based on Michal s statement about discarding the imaginary bins and just using the real ones. While my example might not work for his FFT s
                    Message 9 of 11 , Jun 17, 2002
                    • 0 Attachment
                      I made my comments based on Michal's statement about discarding the
                      imaginary bins and just using the real ones. While my example might
                      not work for his FFT's output format, the method is still valid:
                      a=sqrt(r*r+i*i)/(N/2). If the compiler or library handles this (I
                      know nothing about Python...), that is cool!

                      --- In buildcheapeeg@y..., Jim Peters <jim@u...> wrote:
                      > larryjanow wrote:
                      > > You can't discard the 'imaginary' half of the FFT output, you
                      need
                      > > it to calculate the actual output amplitude at each frequency.
                      >
                      > His python code looks like this:
                      >
                      > def makeSpectrogram(slice):
                      > assert len(slice)==64, "we want 32 bins, so we need 64 samples"
                      >
                      > res = abs(FFT.real_fft(slice))[:-1] # discard 33rd slot (is
                      this okay?)
                      > res = Numeric.floor(res) # round off to integers
                      >
                      > assert len(res)==32, len(res)
                      > return res
                      >
                      > If real_fft is returning 33 values from a 64-sample input, then it
                      > seems to me that it must be returning complex numbers. If so, then
                      > the abs() function should automatically do the square root you
                      > mention (i.e. sqrt(r*r+i*i)).
                      >
                      > Perhaps this needs checking. However, obviously the Python
                      > FFT.real_fft is not returning the imaginary parts as real numbers
                      in
                      > the places you are expecting it to, so your method won't work in
                      this
                      > particular case, and is based on some other FFT algorithm's output
                      > layout.
                      >
                      > Jim
                      >
                      > --
                      > Jim Peters (_)/=\~/_(_) jim@u...
                      > (_) /=\ ~/_ (_)
                      > Uazú (_) /=\ ~/_ (_)
                      http://
                      > B'ham, UK (_) ____ /=\ ____ ~/_ ____ (_)
                      uazu.net
                    • Michal Wallace
                      ... Larry: I was wondering how I d figure out amplitude. Thanks! :) Here s the deal. The FFT module has several functions avaliable:
                      Message 10 of 11 , Jun 17, 2002
                      • 0 Attachment
                        On Mon, 17 Jun 2002, larryjanow wrote:

                        > --- In buildcheapeeg@y..., Jim Peters <jim@u...> wrote:
                        >
                        > > If real_fft is returning 33 values from a 64-sample input, then it
                        > > seems to me that it must be returning complex numbers. If so, then
                        > > the abs() function should automatically do the square root you
                        > > mention (i.e. sqrt(r*r+i*i)).
                        > >
                        > > Perhaps this needs checking. However, obviously the Python
                        > > FFT.real_fft is not returning the imaginary parts as real numbers
                        > > in the places you are expecting it to, so your method won't work in
                        > > this particular case, and is based on some other FFT algorithm's output
                        > > layout.
                        >
                        > I made my comments based on Michal's statement about discarding the
                        > imaginary bins and just using the real ones. While my example might
                        > not work for his FFT's output format, the method is still valid:
                        > a=sqrt(r*r+i*i)/(N/2). If the compiler or library handles this (I
                        > know nothing about Python...), that is cool!


                        Larry: I was wondering how I'd figure out amplitude. Thanks! :)

                        Here's the deal. The FFT module has several functions avaliable:

                        http://pfdubois.com/numpy/doc/FFT/FFT.py.html

                        I used real_fft because the sound processing in python
                        article said the real and imaginary parts would be mirror
                        images. I didn't check that assertion or try to see if it
                        worked in my case. I just chopped it off and saw that the
                        right thing seemed to be happening. But there's a plain
                        fft() function that returns a list of the same length as the
                        input.

                        I should say that complex numbers are a primitive type in
                        python, so the fft() function returns an array of complex
                        numbers, and you can take the .imaginary attribute or the
                        .real attribute to get an array with only the axis you want.

                        I've got some work-related stuff I need to get done today,
                        but I should have some time to play around with this
                        tomorrow.

                        If anyone wants to tak a look in the meantime, be my guest.

                        Cheers,

                        - Michal http://www.sabren.net/ sabren@...
                        ------------------------------------------------------------
                        Switch to Cornerhost! http://www.cornerhost.com/
                        High Powered Hosting - With a Human Touch. :)
                        ------------------------------------------------------------
                      Your message has been successfully submitted and would be delivered to recipients shortly.