Shop OBEX P1 Docs P2 Docs Learn Events
Heater's Fast Fourier Transform. - Page 3 — Parallax Forums

Heater's Fast Fourier Transform.

13»

Comments

  • Heater.Heater. Posts: 21,230
    edited 2011-08-29 03:23
    Ah, OK, that look's better. I'll try and find time to play with this more later today.

    Is there a reason why you have chosen the Hamming window as opposed to all the other possibilities?
  • Andrey DemenevAndrey Demenev Posts: 377
    edited 2011-08-29 06:17
    No special reason except that I have seen this type of window function used in different applications. I think there is enough room to implement other type of windows.

    By the way - I think I need more testing with that window before that code can go into new version



    Heater, I am observing some strange thing. To be sure that is not me, I reverted to your original code. I sample a 1500 Hz sine wave generated by computer's sound card, with a sigma-delta ADC on Propeller, with 48 kHz sampling rate. Then I apply 256 point FFT. Both sampled data and FFT results are sent to serial port. Look at these graphs showing sampled waveform and FFT. They were made by sampling same input signal, one after another. The waveforms are identical (except phase), and visually look pretty clean. But one graph shows that signal has highest magnitude in 8th FFT bin (48000/256=187.5, 1500/187.5=8), and another has a peak in 24th bin (3rd harmonic) !

    good.png


    bad.png


    And while the first graph looks better, I still cannot believe that third harmonic of this signal is only -12dB relative to first. I imported the data sampled by Propeller into Audacity - and this is what it shows, looks more like truth:

    audacity_spectrum.png


    Note that 3rd harmonic is about 40 dB below first, and also there is 2nd harmonic about same level as 3rd.

    Any thoughts?
    725 x 545 - 42K
    725 x 545 - 41K
    946 x 487 - 33K
  • Mark_TMark_T Posts: 1,981
    edited 2011-08-29 08:03
    Hamming windows aren't the best - Kaiser-Bessel windows are significantly better at producing clean looking spectrum when the frequencies are not exact match to the frequency points of the FFT. The Kaiser-Bessel window is also tunable as there is a parameter you can use to adjust its performance. Well that's according to my copy of "Digital Signal Processing - a Practical Approach" - and borne out by my experiments with coding FFTs. If you are only dealing with 8 bits it probably doesn't matter though ;)
  • Heater.Heater. Posts: 21,230
    edited 2011-08-30 02:19
    Andrey,

    Not sure what is up there.

    It is possible something is overflowing. Are your plots showing actual amplitudes? You have a big DC bias on there. What happens when you reduce the signals and/or remove the bias?

    How did you get to make a 256 point FFT, what happens when you use 1024?
  • Andrey DemenevAndrey Demenev Posts: 377
    edited 2011-08-30 04:00
    I do not see reasons why FFT should not work with 256 points. I could expect performance degradation, but not that bad, to totaly unusable results.

    There was very small DC bias - the waveform Y axis is on the right.

    Attached is a test that illustrates the problem. Set FFT size to 256 or 512, and play with PHASE constant in 0 - $800 range - you will see FFT result changing from something resembling the truth to totally wrong. With 1024 points, results seem to be correct regardless of PHASE
  • Heater.Heater. Posts: 21,230
    edited 2011-08-30 05:12
    OK. Should not be overflowing then.

    If it works correctly with 1024 point but fails on 256, 512 my immediate suspicion is that something is going wrong with the "twiddle factors". That is to say get_cos and get_sin are going wrong and not understanding the different ranges required. I'm a bit tied up at the moment so no time to look at it. Perhaps later today.
  • Andrey DemenevAndrey Demenev Posts: 377
    edited 2011-08-30 05:17
    Good luck with it!

    My printer had recently been busy, and now I have few papers about fast Hartley transforms. Fascinating reading before going to bed :)
  • Heater.Heater. Posts: 21,230
    edited 2011-08-30 05:25
    I recently tried to make a C version of the FHT written in RATFOR in Ullman's paper. Never did get it to work.
  • Andrey DemenevAndrey Demenev Posts: 377
    edited 2011-08-30 05:49
    Same here, I failed to make some working code based on Ullman's paper. I am now studying Bracewell's patent
  • HemPeterHemPeter Posts: 28
    edited 2012-01-04 08:20
    Hi Guys,

    I just wanted to say that this thread has been of enormous help to me! I'm currently trying to write an audio tuner with the propeller for my dissertation. I've been spending the last few months getting my head around FFTs, windowing, and all the other goodies involved with this project.
    I've managed to get a sampler to take mic input and send between 9-14 bit samples to between 1 and 4 FFTs to comfortably chug away and display plots on the screen, thanks to Beau's example FFT on the wiki.
    I've been so impressed with the efforts of everyone on this forum! I love the sense of openness of the propeller community. I first started with a propeller at the beginning of 2010 when I was on exchange in Florida. I brought the demo board back to the UK with me and told my Professor about it. He suggested I use it as a final year project and this is the application that I chose.
    I only just discovered the FHT and Constant-Q transforms yesterday, and I fear that I've been heading in the wrong direction for audio. My goal is to be able to tune a Kalimba (finger piano) by playing a tine and displaying the fundamental note and the percentage offset. The problem is the lower frequencies. The range of the Kalimba is B3 to D6.
    So the difference between A3 and A#3 is 13.08Hz, between A#3 and B3 is 13.86Hz, and between B3 and C3 is 14.68Hz. Whereas in the higher range the difference is in the region of 62.23Hz and 65.93Hz
    Note         Frequency
    A#3/Bb3       233.08
    B3            246.94
    C4            261.63
    ...
    C6            1046.50
    C#6/Db6       1108.73
    D6            1174.66
    
    I want decent frequency resolution (at least 1Hz, but preferably ~0.5) of the lower frequencies, whilst being able to detect the highest (only ~1.2KHz) with a 1024 point FFT. I'm a little unsure of the sample rate equation. I gather that I would need a sample rate of double the highest frequency, so 2.4KHz. But then if it's 2400/1024 then the frequency bin resolution would be 2.34Hz. This isn't really enough to be able to discriminate more than ~5.5 steps between notes.. and would that be enough to get the true note for B3?

    Then there's the grey area in my knowledge about band pass filters and envelopes to eliminate the unwanted frequencies and noise.

    Any help or pointers you guys can give you be much appreciated. From advice on windowing or to whether Constant-Q or FHT is the way to go? I can put up the working code I have already if anyone's interested. It doesn't use heater's fft yet, but if you guys think its the way to go it wouldn't be hard to slot it in to the existing project.

    Cheers all
    Pete
  • Dave HeinDave Hein Posts: 6,347
    edited 2012-01-04 09:30
    The number of frequency bins for an FFT is half the tranform size, so a 1024 FFT has 512 frequency bins. If you sample at 2400 Hz the highest frequency you can represent (in theory) is 1200 Hz, so each frequency bin is 1200/512 = 2.34 Hz, which is the same number you computed. If you really want to analyze a tone at 1200 Hz you need to sample at more than twice that frequency. If you're only concerned about the fundamental frequency compoinent I would recommend sampling at 3 or 4 times that frequency. If you want to include the upper harmonics you will need to sample at more than twice the fequency of the upper harmonic you want to include.

    Keep in mind that you need to lowpass filter the signal to less than half the sampling rate to prevent the harmonics above the Nyquist frequency from aliasing to lower frequencies.

    So if you sample at 3 times 1200 Hz and you need 0.5 Hz resolution, you would need a Fourier transform with more than 3600/.5 = 7200 points. This would require an 8192 point FFT.

    EDIT: You could get finer resolution with a smaller FFT if you can determine that the signal is at a low frequency, and subsample the signal before doing the FFT. Of course, there are many other ways to do this without a transform, such as implementing a PLL or implementing a frequency counter.
  • HemPeterHemPeter Posts: 28
    edited 2012-01-04 12:30
    Hi Dave, Thanks for getting back to me and confirming my calculations.
    I was thinking, since one can parallelise the FFT's, would it be silly to take a sample window at a low rate, feed it to 1 FFT, then sample at a faster rate, and let another one handle that? To get the overall gauge of the fundamental, then focus in on one rate once a frequency above a threshold had been detected? Would it create all sorts of problems by changing the sample rate dynamically?
    I'm intrigued by the PLL and frequency counter suggestions. Would it need to be a very pure note for that to work? The resonance of the kalimba is quite temperamental. Some times you get a real twangy sound, like a ruler on a desk. Other times a nice clean note.
  • Dave HeinDave Hein Posts: 6,347
    edited 2012-01-04 13:40
    You may be able to determine the frequency to a higher precision than Fs/N (Sample frequency/FFT size) by comparing the amplitude of the frequency coefficients around the peak value. An isolated peak value would represent a pure sine wave at that frequency. However, if the peak values were spread over two adjacent frequency coefficients the actual frequency would be somewhere in between the two coefficient frequencies.

    The other possibility is to resample the signal as I mentioned above. You would have to apply a lowpass filter before downsampling.

    The PLL or frequency counter would probably require that determine the approximate frequency first. This would be used to set a center frequency for the PLL. The frequency counter would need a lowpass or bandpass filter to remove the harmonics based on the approximate frequency value.
  • HemPeterHemPeter Posts: 28
    edited 2012-01-04 14:30
    OK, that makes sense. So you can use the difference in the relative amplitudes of adjacent bins to determine the given frequency. So if there was a non integer frequency that I was aiming to detect I'd have to look for a certain ratio in 2 particular bins?
    Here's my stab at what process would be:
    Get 1024 audio samples.
    Perform a low pass filter.
    Apply window.
    Run FFT.
    Look for highest peak.
    Check note look up table for expected adjacent peak ratios.
    Determine percentage offset between current note and the next note.
    Display on the screen.

    Is that the right way round for LP filter and windowing?

    I'm not entirely sure what process you mean what you talk about subsampling and down sampling..
  • Dave HeinDave Hein Posts: 6,347
    edited 2012-01-04 14:44
    The low-pass filter is done in the analog domain before you sample the signal. You can probably get by with a simple RC filter with a cutoff around 1 KHz.
  • nglordinglordi Posts: 114
    edited 2012-01-04 14:54
    Andrey, Heater, HemPeter:

    There are Spin & PASM FHT versions based on Ullman's paper posted on the forum.
    Search for Fast Hartley Transform.

    NickL
  • Heater.Heater. Posts: 21,230
    edited 2012-01-05 04:16
    nglordi,

    Yes I have seen the FHT. Haven't had the time to play with it yet. When I tried to translate Ullman's RATFOR FHT to C a while back I could not get it to work. Did you find any issues with Ullman's RATFOR code?

    HemPeter,

    Interesting project, never knew a finger piano was called a Kalimba. How does one tune a Kalimba? Do you have to file down the keys?

    I think that PASM FFT on the wiki is by Ale not Beau. The header contains refers to "Pacito" as the author which is in Ale's sig.

    The Heater FFT came about because after many years wondering about it I finally understood (with help from the great folks on this forum) how the FFT works and could write one from scratch. Having done that I C of course I had to do it in PASM.

    I have never played with the FFT on wiki so I don't know how it compares in accuracy and or performance, I'm sure it's fine. As it stand that FFT is licensed under the GPL where as the Heater FFT is MIT licensed.

    Not sure I have any suggestions for your project yet. Seems to me that if you want to resolve frequency down to 0.5Hz you need to have two seconds worth of samples. But then as pointed out perhaps it's enough to look at the peaks spanning a few frequency bins and interpolate a bit.
  • nglordinglordi Posts: 114
    edited 2012-01-07 18:52
    Heater:

    With respect to Ullmann's code, there was a problem with the bit reversal loop. I made the following changes
    in the code.

    butdis = 1 (not 0)
    if (k =< butdis) (not <=)
    while (......... && butloc >1 (not >=)


    NickL
  • HemPeterHemPeter Posts: 28
    edited 2012-01-08 21:00
    Dave,
    I found the FIR2PASM: Automatic FIR Filter Code Generator thread. Do you think that would be adequate for implementing the low pass filter? I would like to ideally like to use the mic on the demo board.
    I have a vague plan of how to plug everything together (control, LP filter & windowing sampler, vga, serial debug, fft/fht x 1-4)
    Forgive me if these sound like noob questions. I'm still relatively new to all this.
    With a low pass filter, windowing and sampling at 500Hz with a 1024 point FFT, would there be any reason that I could not discriminate a pure A3 (220Hz) at 0.5Hz resolution? Would the very low level noise require a band pass filter? Would there not be a single spike in a single bin if all that was present was the pure note?

    nglordi,
    Thanks for the link, would you say FHT would be the way to go with this? I haven't looked into it yet.

    Heater,
    Tuning a Kalimba is relatively difficult. It consists of 17 different length flat pieces of metal (tines) held in place at the top by a bar that's screwed down. The easiest method I've found for tuning is pliers. It requires some to-ing and fro-ing to get the tine to move at all, and then any further movement causes it to go out of tune fairly easily.
    Cheers for the wiki FFT credit correction. Did you have to register the Heater FFT code before putting an MIT licence on it? I commend you for managing to get your head round it! Just the terminology has been a brain twister for me.

    I'm wondering if I should battle through with the challenge of understanding the FFT at this point. If I could achieve higher accuracy with a dynamic frequency counter or the FHT or Constant-Q transform I could explore those options. I have to have a prototype ready in some way, shape or form by the 25th of Jan so pressure's beginning to mount a bit.

    If you guys know if any better thread this topic of discussion should be moved to let me know, or I suppose I could start me own? =)
    Ahh, coding in the morning ;)
  • Heater.Heater. Posts: 21,230
    edited 2012-01-09 04:15
    HemPeter

    Re: Filtering.

    The name of the game here is to keep frequencies above half of your digital
    sampling rate out of the digital system else things start to go badly wrong.
    This applies to the FFT and to the FIR filter and any other digital domain
    stuff.

    This means that the filtering should exist in the analog domain prior to being
    sampled.

    The reason for all this is that frequencies above half of your sampling rate
    will appear in the samples as a signal at some lower "alias" frequency. Google
    around for "aliasing" and "Nyquist–Shannon sampling theorem".

    Having said that if, for example, you were doing this on your PC with a 44KHz
    sampling rate, and max input frequency of 20KHz. All OK so far. Then for some
    reason you want your FFT to only sample 1024 times per second I guess you could
    use a FIR to remove high frequencies from your samples prior to hitting the
    FFT.

    Now as for that 220Hz. Lets sample the analog at 1024 samples per second. Our
    maximum allowed analog input frequency is 512Hz. Run that into an FFT or FHT or
    1024 samples. The FFT/FHT output bins are now 1Hz to 512Hz. With a 1Hz
    resolution.

    Seems to me to get to 0.5 Hz resolution you need to sample over two seconds. So
    lets sample at 512Hz. Our max allowed input frequency is now 256Hz. Take 1024
    samples into FFT/FHT and we have output frequency bins of 0.5Hz to 256Hz. That
    220Hz just about fits.

    Ahhh, this is giving me headache...

    No you don't to register anything to slap whatever license terms you like on
    your creations. Whilst those terms are called the MIT license, as that is where
    they originated, when I or you use them there is no connection to MIT. Note
    that all code donated to the Parallax object exchange requires licensing that
    way.

    Sure you should try and understand the FFT/FHT but at this point I would not
    try to hard. First you will have to understand the regular Fourier Transform
    mathematics, then the Discrete Fourier Transform (DFT). Concentrate on what the
    input and output means and how it can be used. Whilst we studied the Fourier
    Transform as students it was not until 20 years after first discovering the FFT
    algorithm that I understood how it works!!

    For your prototype I would put together the sampling, FFT/FHT, VGA, serial and
    show that you can identify frequencies in the input. Perhaps use a signal
    generator to feed in clean sine waves to avoid "alias" problems initially. Show
    what happens when using real Kalimba input and be ready to explain how you are
    going to attempt to fix any issues. That's already quite a lot.

    I thing starting your own "kalimba tuning" thread is a great idea. I want to see
    how this pans out.
  • Heater.Heater. Posts: 21,230
    edited 2012-01-09 04:22
    P.S. This all reminds be of watching a couple of guys in my final year undergraduate physics group in the lab trying to measure the frequency of a big brass bell with a digital frequency counter, they were wondering why they got pretty much random numbers out of it. I was wondering what they had learned in 3 years!!!

    You are already many steps ahead...
  • cessnapilotcessnapilot Posts: 182
    edited 2012-01-11 11:11
    HemPeter:

    What you need to achieve that high frequency resolution with small FFT size has been already invented. Look for Zoom FFT on the net. That simple method, however, uses complex numbers and floating point FFT for high enough dynamic range to avoid spurious effects with integers and fixed points. In Zoom FFT you can select a small spectral range to span with the same number of frequency bins as with the original (from zero to max. freq.) However, when you inrease resolution, say 5 times, you have to oversample the signal 5 times before digital filtering. In other words, you cannot overcame Heisemberg's uncertainty principle for time and energy (here frequency) with mathematical tricks.

    The FPU32_FFT driver on OBEX can be easily modified to perform Zoom FFT with the coprocessor. The next version of the driver will certainly contain that option.

    -Istvan
  • bryenbryen Posts: 1
    edited 2012-01-12 14:35
    Folks on this thread may be interested in a DFT algorithm put forward by a former colleague of mine :

    http://www.icita.org/icita2004/abstracts/128-5.htm

    [h=1]128-5:A Cyclic Group Based Prime Number DFT Method for the Calculation of the Cepstrum[/h] Dr. David Tien
    Charles Sturt University, Australia

    Abstract: Since the introduction of the concept of cepstrum by Bogert [1] in the early 60s, the theoretical and computational aspects of cepstrum analysis have been widely used with success in a wide range of signal processing areas. One of the major methods to compute the cepstrum is to employ the technique of Discrete Fourier Transformation (DFT). However, the direct evaluation of an N-point DFT needs about 4N2 real multiplications as well as additions which slows the calculation substantially. Even the Fast Fourier Transformation (FFT), introduced by Cooley and Tukey [2], still requires 2Nlog_2(N) operations. Thus, for reasonably large sample values of N, such existing methods require an inordinate number of computations. This in turn has made the FFT processing slow and the calculation of the cepstrum expensive. To combat the above mentioned problems, in this article we have developed a fast algorithm which will substantially reduce the complexity of the current available cepstrum analysis algorithms.

    Index Terms: Cyclic Group, Prime Number, DFT Method, Cepstrum.


    The PDF is available from the link above, ****but I am no expert (or even beginner) on this Prime Num DFT algorithm, and have never tried to implement or understand it, so I have no enlightening comments to make about it. He claims a processing reduction of up to about 95%... your mileage may vary :)
  • dbpagedbpage Posts: 217
    edited 2012-02-20 13:17
    If you find that bryen's link is broken, the article can be found at: http://researchoutput.csu.edu.au/R/?func=dbin-jump-full&object_id=6815&local_base=GEN01-CSU01
  • HemPeterHemPeter Posts: 28
    edited 2012-06-19 11:32
Sign In or Register to comment.