Help needed with Fast Fourier Transform

[Unscruffed]

A very generalized title, but quite frankly I'm lost altogether with this.

Here's what I want to do: Monitor the soundcard input, analyze it, and get the frequency of the input signal in Hertz (as a Double).

All the reading I've done seems to point to using FFT to acomplish this. However all the code samples I find, seem to be based around spectrum analyzers or displaying data in other similar ways. I can't find any examples of how to simply get the frequency of the signal.

The code I have for monitoring the soundcard input returns a buffer via an event (as an array of Integers) every 100 milliseconds or so. But I have no idea of how to analize this buffer.

I've looked at tons of code examples from PSC and elsewhere, but they are all far more complex than what I need, and are so poorly documented and/or formatted to the extent that I can't even follow the flow of the code to be able to learn anything at all from them.

Can anyone help out with a basic example of how to achieve what I'm trying to do?

The buffer is 1 channel @ 22050 sample rate, and is provided via an event as follows:

Private Sub WaveMonitor_GotData(Buffer() As Integer, Length As Long)
    ' What to do here?
End Sub

Thanks in advance for any help.

Edit:
FYI, the monitoring is done using a WaveRecorder class (which I got from an example on PSC from memory).
It uses the waveIn... functions in WINMM.dll. DirectX is not used at all.

 

[Unscruffed]

This is the simplest code that I've been able to find so far, and for what does, it does well. But just like all the other so-called "frequency detection" examples I've found, this code relies on comparison, not detection - which are two entirely different things.

Public Function Power(Wave() As Single, Frequency As Double, SampleRate As Long) As Single
    ' Wave() is filled with samples with the range -1 to +1
    ' Frequency is the frequency to measure, in cycles per second.
    ' SampleRate is the audio recording rate - usualy 44100 or 48000
    ' Uses discrete fourier technique to measure power amplitude of specific
    ' frequencies in audio data.
    Dim lCount As Long
    Dim sMulti As Single
    Dim sReal As Single
    Dim sImag As Single
    Dim i As Long
    Const Pi = 3.14159265
    lCount = (1 + UBound(Wave))
    sMulti = (Frequency * 2 * Pi / SampleRate)
    For i = 0 To (lCount - 1)
        sReal = (sReal + Wave(i) * Cos(i * sMulti))
        sImag = (sImag + Wave(i) * Sin(i * sMulti))
    Next
    sReal = (sReal / lCount)
    sImag = (sImag / lCount)
    Power = (sReal * sReal + sImag * sImag)
End Function

It's simply not feasible to test for every frequency in the audible spectrum. If you did that, and you required accuracy of 1Hz, then you would need to call this function around 20,000 times every time the buffer presents itself.

Comparison is simply not an option here - I need proper detection, but can't find one single example of how to do it.

Heaven doesn't want me, and Hell's afraid I'll take over!

 

[Hypetia]

This will very much depend upon the shape and specification of your input signal. If your signal is monotone, and has fixed amplitude, then it will be very easy to detect. However, if your signal is more complex, for instance, contains multiple frequencies of varying amplitude over time, than you might really need some kind of Fourier transform to analyze the data. From your post, I get the feeling that your signal is monotone.

In your WaveMonitor event, you need to write a unidirectional zero-crossing detector to count the number of troughs or crests of the incoming signal sample. Once you get the count, you can compute the frequency as follows.

Assuming the 'Length' parameter represents the number of samples in Buffer().

Sample Rate = 22050
Sample Length (in seconds) = Length / Sample Rate (Type: Double)
Frequency = Count / Sample Length

Assuming the event occurs every 100ms, the Length should be typically 2205 if all samples are returned from the last event. (Sample Rate * 0.1s = 2205).

If the zero crossing count is, say 500, then frequency will be calculated as follows.

Sample Rate = 22050
Sample Length (in seconds) = 2205 / 22050 = 0.1
Frequency = 500 / 0.1 = 5kHz

You should know exactly what does Length means and write your code accordingly.

You may also need to implement a hysteresis filter instead of simple zero-crossing detector to filter out noise or any weak frequency from the signal. That will require knowledge of how Buffer() is filled and interpreted. Check the documentation of WaveRecorder class.

 

[Unscruffed]

Thanks for the reply. I apologise for my oversight of not providing enough info about the source signal - that obviously would have helped.

What you've suggested, is essentially what I'm doing now. I've got this working for one purpose, but not another. So I'll try to explain the exactly what I'm trying to do.

The program serves two purposes:
1. a guitar tuner, and
2. an instrument for measuring the integrity of the signal (frequency drift).

Obviously, I have no idea how much anyone knows about guitars, so I'll try to give as much detail as I can.

For tuning, there's two different methods that are commonly used to generate the signal:
1. playing the open string (not fretting it), or,
2. playing the harmonic at the 5th fret (by gently touching the string above the 5th fret and then releasing it - without actually fretting it).

These two different approaches, result in very different signals. Playing any open or fretted note, results in a signal which contains the note frequency plus harmonics. The actual frequency of the target note is the most dominant frequency in the signal. You therefore need to detect the dominant frequency.

A harmonic at the 5th fret however, essentially produces a perfect sine wave which is 2 octaves above the tuned frequency of the string. With a bass guitar using standard tuning for example, the open string notes are E1, A1, D2 and G2. The harmonics at the 5th fret produce sine waves at E3, A3, D4 and G4 respectively. So detecting the frequency in this scenario, is as simple as measuring the number of samples between zero crossings which occur in the same direction (which is every second crossing), and doing some simple calculations. (Exactly as you described in your reply.)

The integrity check however, can only be performed by playing the open note, letting it ring, and then measuring how far the frequency drifts. Playing an open note causes the most vibration, and therefore produces the highest amount of tension on the strings and the neck. Being able to measure this is useful for two reasons...

Firstly, the older or cheaper the stings are, the more drift you will get because of stretching. Admittedly, this is not overly useful in general terms as a player, because you can hear when the strings are on the way out anyway. But it can be useful for checking the performance of the strings, or comparing things such as different brands and gauges of strings.

The second and more important purpose, is that with an older or cheaper guitar, or a guitar that has been stored for a long period of time with the strings at normal (tuned) tension, there can also be movement in the neck that causes some of this drift. Measuring the amount of drift therefore, particularly with bass guitars which use much thicker and stronger strings, can be useful for checking the integrity of the neck itself. A guitar neck is slightly bowed when under tension, and it straightens as you loosen the strings and/or remove them. This bow, is controlled by what's called a "truss rod", which is a metal rod inside the neck itself, and runs it's full length. Excess frequency drift could therefore indicate that the truss rod needs adjusting, or even more serious issues, such as a truss rod or neck that should actually be replaced. (High-end guitar manufacturers actually perform these types of tests during R & D to evaluate neck designs.)

So to sum it all up, for tuning using harmonics, I can use the simple zero crossing method to measure the frequency, and this is very accurate because the signal is a near perfect sine wave.

But for tuning by playing the open strings, or for performing the integrity check, I need to detect the dominant frequency in a non-perfect signal. Furthermore, for the integrity check, I also need to measure the frequency drift over time as the signal decays - which will be straight forward once the detection problem is solved.

So, all my research so far has led me to Fast Fourier Transform (FFT). From what I understand, this converts data from a "time domain" to a "frequency domain". But what I can't find, is how to use this transformed data to determine the dominant frequency. In the example code I posted above, the routine returns the power of a known frequency. All the examples I've found so far do similar things - claiming to be "detection" routines, when in fact they are nothing more than "comparison" routines because they rely on a known frequency to test against.

I simply don't understand anywhere near enough about FFT to be able to figure out how to use it - mathematically, it's way over my head. I'm not even sure if FFT is what I actually need - since every search for "frequency detection" points in that direction, but I can't find any example of "real" detection taking place.

 

[strongm]

Hmm ...

ALGLIB contains FFT methods. There's free VB6-compatible source code for this maths library available here.

Whether it is suitable for your purpose - or actually quick enough - I can't say.

 

[Unscruffed]

G'day Strongm, been a long time. Hope all is well.

I already have that zip. It contains a whole bunch of *.bas modules that I couldn't get to work. They heavily cross reference each other (with types and so on), and I had nothing but problems when I tried to use it. They don't even use Option Explicit - which is very poor practice IMO. I couldn't really learn anything from them.

I have however tracked down one FFT module which (after some formatting) is clean and seems to work, but I'm questioning the way that it calculates the frequency.

There's two main routines. The first is the actual FFT function:

' FourierTransform '
Public Sub FourierTransform(NumSamples As Long, RealIn() As Double, ImageIn() As Double, RealOut() As Double, ImagOut() As Double, Optional InverseTransform As Boolean = False)
    Dim i As Long
    Dim j As Long
    Dim k As Long
    Dim n As Long
    Dim AngleNumerator As Double
    Dim NumBits As Byte
    Dim BlockSize As Long
    Dim BlockEnd As Long
    Dim DeltaAngle As Double
    Dim DeltaAr As Double
    Dim Alpha As Double
    Dim Beta As Double
    Dim TR As Double
    Dim TI As Double
    Dim AR As Double
    Dim AI As Double
    If InverseTransform Then
        AngleNumerator = -2# * Pi
    Else
        AngleNumerator = 2# * Pi
    End If
    If (IsPowerOfTwo(NumSamples) = False) Or (NumSamples < 2) Then
        Call MsgBox("NumSamples is not a positive integer power of two.", , "Error!")
        Exit Sub
    End If
    NumBits = NumberOfBitsNeeded(NumSamples)
    For i = 0 To (NumSamples - 1)
        j = ReverseBits(i, NumBits)
        RealOut(j) = RealIn(i)
        ImagOut(j) = ImageIn(i)
    Next
    BlockEnd = 1
    BlockSize = 2
    Do While BlockSize <= NumSamples
        DeltaAngle = AngleNumerator / BlockSize
        Alpha = Sin(0.5 * DeltaAngle)
        Alpha = 2# * Alpha * Alpha
        Beta = Sin(DeltaAngle)
        i = 0
        Do While i < NumSamples
            AR = 1#
            AI = 0#
            j = i
            For n = 0 To BlockEnd - 1
                k = j + BlockEnd
                TR = AR * RealOut(k) - AI * ImagOut(k)
                TI = AI * RealOut(k) + AR * ImagOut(k)
                RealOut(k) = RealOut(j) - TR
                ImagOut(k) = ImagOut(j) - TI
                RealOut(j) = RealOut(j) + TR
                ImagOut(j) = ImagOut(j) + TI
                DeltaAr = Alpha * AR + Beta * AI
                AI = AI - (Alpha * AI - Beta * AR)
                AR = AR - DeltaAr
                j = j + 1
            Next
            i = i + BlockSize
        Loop
        BlockEnd = BlockSize
        BlockSize = BlockSize * 2
    Loop
    If InverseTransform Then
        'Normalize the resulting time samples...
        For i = 0 To NumSamples - 1
            RealOut(i) = RealOut(i) / NumSamples
            ImagOut(i) = ImagOut(i) / NumSamples
        Next
    End If
End Sub

This one seems straight forward to use:
[ol 1]
[li]For RealIn(), convert my sample data from an integer array to a double array, with a size that is a power of 2.
(100ms @ 22050 = 2205 samples available each time my event fires. I therefore need to use 2048 samples.)[/li]
[li]I don't have ImageIn() data, so according to the comments in the code, I should pass in an array of 2048 zeros.[/li]
[li]Initialize two more arrays of the same size (2048) to accept the output.[/li]
[/ol]

After calling the function, as far as I can tell, the next step is to search the RealOut() array, and get the index of the highest value. This is then passed to the next routine - which is the one I have a serious problem with:

' FrequencyOfIndex '
Public Function FrequencyOfIndex(NumberOfSamples As Long, ByVal Index As Long) As Double
    If Index >= NumberOfSamples Then
        FrequencyOfIndex = 0#
        Exit Function
    ElseIf Index <= NumberOfSamples / 2 Then
        FrequencyOfIndex = CDbl(Index) / CDbl(NumberOfSamples)
        Exit Function
    Else
        FrequencyOfIndex = -CDbl(NumberOfSamples - Index) / CDbl(NumberOfSamples)
        Exit Function
    End If
End Function

There's no need to even look at this code - just look closely at the parameters: NumberOfSamples will always be 2048, and Index will always be between 0 and 2047 inclusive.

This means that this function can only return one of 2048 values - which can't possibly in any way be accurate enough.

As far as I can gather, this could only provide a "frequency band", which has a bandwidth that is directly proportional to the number of samples. Considering an audible bandwith of 20 to 20,000 Hz, this function can only provide accuracy of ((20,000 - 20) / 2048) - which gives an absolute minimum bandwith of approx 9.75 Hz. (I already have a guitar tuner program which has an accuracy of 0.001 Hz, but unfortunately the source is not publicly available, and it was probably written in C anyway.)

If this is the best that FFT can do when it comes to frequency detection, then it's totally inadequate. I need something else entirely.


Heaven doesn't want me, and Hell's afraid I'll take over!

 

[Hypetia]

Two things.

1. You said you already tried the zero-crossing technique. What problems did you face with that method? If the results are not accurate enough (like getting a frequency higher than actual in result), you should implement the hysteresis filter. It would eliminate noise and low amplitude harmonics.

2. If you go for FFT, the results will always be discrete as FFT is derived from DFT. This means you will never get the resulting amplitude as a continuous function of frequency. The FFT will return amplitude at discretely separate frequencies and spacing between those frequencies will depend upon the sample rate of the input signal, which is 22kHz in your case.

The math involved in frequency detection is simple. Since your sample rate is 22kHz, you'll not be able to detect frequencies above 11kHz, whether you use FFT or not. Remember Nyquist theorem? In order to cover whole audio range, you should sample the data at least at 40kHz. The closest standard sample rate is 44.1 kHz --- twice the sample rate you are currently using.

 

[Unscruffed]

@Hypetia

Re #1: That's already been explained in detail, however filtering out the harmonics might be the solution. I'll look into that.

Re #2: That's pretty much what I've discovered. FFT tells us the strength of a series of frequency bands in the signal. The number of bands and their bandwidth is directly related to the number of samples that the calculation is based upon. In this case, I need a much narrower bandwith than FFT can provide.

Re SampleRate: The highest frequency that I need to detect would be the 5th Fret Harmonic for the 6th (thinest) string on a standard guitar. This would be note E6 (1318.51 Hz). The sample rate of 22KHz is therefore fine. In fact, I could go down to 11KHz and still have plenty of room to spare.

Below are the frequency tables for guitars using standard tuning. (I've used a 5 string bass here because of the additional lower B string.)

[pre]
Bass Guitar Tuning Frequencies
Open Note 5th Fret Harmonic
Note Freq (Hz) Note Freq (Hz)
B0 30.868 B2 123.471
E1 41.203 E3 164.814
A1 55.000 A3 220.000
D2 73.416 D4 293.665
G2 97.999 G4 391.995

Standard Guitar Tuning Frequencies
Open Note 5th Fret Harmonic
Note Freq (Hz) Note Freq (Hz)
E2 82.407 E4 329.628
A2 110.000 A4 440.000
D3 146.832 D5 587.330
G3 195.998 G5 783.991
B3 246.942 B5 987.767
E4 329.628 E6 1318.510
[/pre]


Heaven doesn't want me, and Hell's afraid I'll take over!

 

[strongm]

This looks like it might be promising:

https://www.planet-source-code.com/vb/scripts/ShowCode.asp?txtCodeId=64918&lngWId=1

 

[Unscruffed]

I've already checked all the examples at PSC that I can find. Nothing I've found can return the dominant frequency with any accuracy. The closest I've found is the code I posted earlier.

I'm currently following up Hypetia's idea of filtering out the harmonics. So now I'm having trouble finding filter examples in VB.

Heaven doesn't want me, and Hell's afraid I'll take over!

 

[Unscruffed]

Okay. I found a Butterworth filter here.

I've done some experimenting, and this filter works very well. So the next step is to be able to calculate on the fly the cutoff frequency to use for the low pass filter. Here's what I'm thinking might work:

[ol 1]
[li]Perform FFT on the buffer.[/li]
[li]Get the 2 highest frequencies from the FFT results. (The lower of these 2 frequencies should be close to the dominant frequency, and the higher one should be close to the first harmonic.)[/li]
[li]Calculate the frequencey that is half way between the 2 frequencies from step 2.[/li]
[li]Apply a low pass filter using the frequency from step 3 as the cutoff.[/li]
[li]Detect the dominant frequency in the modified buffer by using zero crossing or peak methods.[/li]
[/ol]

You guys think this might work?

Obviously, I only need to apply this process if the source is a played note. (The 5th fret harmonics produce a clean enough sine wave without the need for any of this.)

Heaven doesn't want me, and Hell's afraid I'll take over!

 

[Hypetia]

If you get the dominant frequency in step 2, why do you proceed further? I assume you do so because the frequency returned by FFT is not accurate enough. Is that right?

I did not mean to implement a full-fledged DSP filter. I just wanted you to add a hysteresis window to your zero-crossing code. See the code below which explains the idea.
___
[pre]

Dim Wave(100000) As Single
Const SampleRate = 22050
Private Sub Form_Load()
AddSignal 500, 1 'main signal (500Hz, Amplitude = 1)

Debug.Print Frequency(0) 'calculate frequency without hysteresis. Pure signal, no noise. Result correct.

AddSignal 1234, 0.5 'add some noise (1234 Hz, Amplitude = 0.5)

Debug.Print Frequency(0) 'calculate frequency without hysteresis, result incorrect due to noise.
Debug.Print Frequency(0.5) 'calculate frequency with hysteresis, result corrected.

End
End Sub

Private Sub AddSignal(Frequency As Double, Amplitude As Double)
Dim N As Long
Const Pi = 3.14159265358979
For N = 0 To UBound(Wave)
Wave(N) = Wave(N) + Amplitude * Sin(2 * Pi * Frequency * N / SampleRate)
Next
End Sub

Private Function Frequency(Hysteresis As Double) As Double
Dim N As Long, Flag As Integer, Count As Long
Flag = 1
For N = 0 To UBound(Wave)
If Wave(N) * Flag > Hysteresis Then
Count = Count + 1
Flag = -Flag
End If
Next
Frequency = Count * SampleRate / (UBound(Wave) + 1) / 2
End Function[/pre]
___

This is just a demo code, and results are pretty clean as the waveform is not a real-world signal, but generated from code. You might be able to test actual performance when you try it on your input signal. You might also want to tune the hysteresis window size (0.5) to get suitable results.

>You guys think this might work?
I didn't check the code myself, so can't say anything for sure.