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

need help with DSP

Expand Messages
  • Michal Wallace
    Hey all, Well, I ve got some code that generates some fake data and runs the FFT on it to produce a spectogram, but I need some help. The good news is that
    Message 1 of 11 , Jun 15, 2002
    • 0 Attachment
      Hey all,

      Well, I've got some code that generates some fake data and
      runs the FFT on it to produce a spectogram, but I need some
      help.

      The good news is that even with the FFT it still runs quite
      quickly. The bad news is that it doesn't show what it ought
      to be showing.


      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. The code is
      short, so the bugs can't be in too many places... :)

      Can someone take a look? The new code is here:

      http://cvs.sabren.com/sixthdev/cvsweb.cgi/openeeg/proto.py?rev=1.2


      I used this article on analyzing sound files with
      python as my gude, so it might help:

      http://www.onlamp.com/pub/a/python/2001/01/31/numerically.html

      The dislin package uses for plotting is free, cross-platform
      and very quick, btw... Maybe it's possible to use it for an
      oscilliscope display? Or even a better spectrogram. By default,
      though, it pops up its own window. Maybe there's some way to
      just generate the graph to a sprite and blit it onto the screen
      ourselves as we animate?


      Anyway, if someone with DSP skill can help me figure out
      what's wrong with my code, that would really be great!
      [Even if you don't know python, the code should be pretty
      straightforward]


      Thanks,

      - Michal http://www.sabren.net/ sabren@...
      ------------------------------------------------------------
      Switch to Cornerhost! http://www.cornerhost.com/
      High Powered Hosting - With a Human Touch. :)
      ------------------------------------------------------------
    • 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 2 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 3 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 4 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 5 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 6 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 7 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 8 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 9 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 10 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 11 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.