[vDSP][Signal Processing] Audio Analysis 5 FFT (English ver.)

[vDSP][Signal Processing] Audio Analysis 5 FFT

This article is also readable in
 http://sourcecodedust.blogspot.ch/2015/06/vdspsignal-processing-audio-analysis.html

Here I introduce how to implement Fast Fourier Transform by using vDSP (Acceleration.framework) with C code.


## Process

1. Define FFT setup
2. Extract audio data of FFT frame size
3. Windowing the audio data obtained in the process 2
4. Convert audio data to DSPSplitComplex types.
5. Execute FFT
6. Do scaling

In the fourth process, we find a data type **DSPSplitComplex** that can store both real and imaginary numbers. The structure of the data type of **DSPSplitComplex** is shown as follows,

<pre class="brush: c++;">
typedef struct DSPSplitComplex {
    float  *realp;
    float  *imagp;
} DSPSplitComplex;
</pre>

**DSPSplitComplex** is used for **float** type of real and imaginary numbers. On the other hand, the **double** data type is treated by **DSPDoubleSplitComplex** .


## Window
There are several vDSP functions for making window functions.
Here I introduce some of the most popular window functions.

hanning 

<pre class="brush: c++;">
extern void vDSP_hann_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
</pre>

hamming

<pre class="brush: c++;">

extern void vDSP_hamm_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
</pre>

blackman

<pre class="brush: c++;">
extern void vDSP_blkman_window(
    float       *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
</pre>

etc.
Thoese of all functions introduced above are for **float** data type and are adopted by adding **D** at the end of function name for **double** data type.
For example, hanning winwodw can be written as below,

<pre class="brush: c++;">
extern void vDSP_hann_windowD(
    double      *__vDSP_C,
    vDSP_Length  __vDSP_N,
    int          __vDSP_Flag)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);

</pre>


Now, I show you the main function which handles FFT analysis.

##main()

<pre class="brush: c++;">
#define PAI 3.14159265358979323846264338 // PAI
#define SAMPLE_RATE 44100 //sampling rate
#define LENGTH 4 // 4 seconds
#define AMPLITUDE 1.0// 16bit
#define FREQUENCY 440 // frequency Hz
#define FFTSIZE 2048

int main(int argc, const char * argv[])
{
    //Libsndfile
    SNDFILE *fp;
    SF_INFO sfinfo;
    float *buffer;
    buffer = malloc(FFTSIZE*sizeof(float));

    memset(&sfinfo,0,sizeof(sfinfo));

    sfinfo.samplerate = SAMPLE_RATE;
    sfinfo.frames = SAMPLE_RATE*FFTSIZE;
    sfinfo.channels = 1;
    sfinfo.format = (SF_FORMAT_WAV | SF_FORMAT_PCM_16);

    if(!(fp = sf_open("Sawtooth.wav", SFM_WRITE, &sfinfo))){
		printf("Error : Could not open output file .\n");

		return 1;

	} 
	
	
    int i;
    for(i=0;i&lt;FFTSIZE;i++){
        buffer[i] = AMPLITUDE * sin((float)FREQUENCY / SAMPLE_RATE * 2 * PAI * i);
    }
    
    

 int ffthalfsize = FFTSIZE/2;
    DSPSplitComplex splitComplex;

    //allocate complex (fft half size)
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);

    float *magnitude = malloc(sizeof(float)*ffthalfsize);

    FFT(buffer,&amp;splitComplex);

    vDSP_zvabs(&amp;splitComplex, 1, magnitude, 1, ffthalfsize);

    float freq_bins = (float)SAMPLE_RATE/FFTSIZE;

    for(i=0;i&lt;ffthalfsize;i++){
        printf(&quot;%0.2f Hz: magnitude = %f\n&quot;,i*freq_bins,magnitude[i]);
    }
            FFT_free(&amp;splitComplex);

    if(sf_write_float(fp, buffer,sfinfo.channels *FFTSIZE)!= sfinfo.channels *FFTSIZE)
        puts(sf_strerror (fp));

    free(buffer);    
    sf_close(fp);

    return 0;

}

</pre>


Here we focus on the place after sin wave is generated.

<pre class="brush: c++;">
    int ffthalfsize = FFTSIZE/2;
    DSPSplitComplex splitComplex;

    //allocate complex (fft half size)
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);  

    float *magnitude = malloc(sizeof(float)*ffthalfsize);

    FFT(buffer,&amp;splitComplex);

    vDSP_zvabs(&amp;splitComplex, 1, magnitude, 1, ffthalfsize);

    float freq_bins = (float)SAMPLE_RATE/FFTSIZE;

    for(i=0;i&lt;ffthalfsize;i++)
        printf(&quot;%0.2f Hz: magnitude = %f\n&quot;,i*freq_bins,magnitude[i]);

    FFT_free(&amp;splitComplex);
</pre>

In the soruce code, the FFT frame size = 2048.
Make Half FFT frame size.

<pre class="brush: c++;">
    int ffthalfsize = FFTSIZE/2;
</pre>

Subsequently, declare **DSPSplitComplex** data and allocate memory spaces for both real and imaginary numbers. The size of both scolars are half of fft frame size.

<pre class="brush: c++;">
    splitComplex.realp = malloc(sizeof(float)*ffthalfsize);
    splitComplex.imagp = malloc(sizeof(float)*ffthalfsize);
</pre>

Next, allocate memory space for storing spectral magnitude. The size of the space is half fft frame size.
<pre class="brush: c++;">
     float *magnitude = malloc(sizeof(float)*ffthalfsize);
</pre>

Finally, we execute FFT. 
I made FFT() user function whose instruction is shown after this sub-section.

**FFT()** takes pointer of **buffer** and **splitComplex** as its arguments.

<pre class="brush: c++;">
FFT(buffer,&splitComplex);
</pre>

When **FFT()** is executed, the result is stored in real and imaginary spaces of **splitComplex**.
We can obtain spectral magnitude by calcurating the power of these numbers.

<pre class="brush: c++;">
vDSP_zvabs(&splitComplex, 1, magnitude, 1, ffthalfsize);
</pre>

**vDSP_zvabs()** calculates sqrt(real * real + img * img) and stores the result in **magnitude**.

<pre class="brush: c++;">
extern void vDSP_zvabs(
    const DSPSplitComplex *__vDSP_A,
    vDSP_Stride            __vDSP_IA,
    float                 *__vDSP_C,
    vDSP_Stride            __vDSP_IC,
    vDSP_Length            __vDSP_N)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);
</pre>

**Stride** is 1.0.
**vDSP_Length** represents the data length of the calculation and so that it should be **ffthalfsize**.

Calculate which FFT bin corresponds what frequency.
<pre class="brush: c++;">
float freq_bins = (float)SAMPLE_RATE/FFTSIZE;
</pre>

Output the result.

<pre class="brush: c++;">
for(i=0;i&lt;ffthalfsize;i++)
        printf(&quot;%0.2f Hz: magnitude = %f\n&quot;,i*freq_bins,magnitude[i]);

</pre>

Here some of the result.

<pre class="brush: c++;">
0.00 Hz: magnitude = 0.000018
21.53 Hz: magnitude = 0.000018
43.07 Hz: magnitude = 0.000019
64.60 Hz: magnitude = 0.000020
86.13 Hz: magnitude = 0.000023
107.67 Hz: magnitude = 0.000026
129.20 Hz: magnitude = 0.000030
150.73 Hz: magnitude = 0.000036
172.27 Hz: magnitude = 0.000044
193.80 Hz: magnitude = 0.000055
215.33 Hz: magnitude = 0.000072
236.87 Hz: magnitude = 0.000096
258.40 Hz: magnitude = 0.000134
279.93 Hz: magnitude = 0.000195
301.46 Hz: magnitude = 0.000301
323.00 Hz: magnitude = 0.000504
344.53 Hz: magnitude = 0.000943
366.06 Hz: magnitude = 0.002103
387.60 Hz: magnitude = 0.006500
409.13 Hz: magnitude = 0.051471
430.66 Hz: magnitude = 0.221126
452.20 Hz: magnitude = 0.202370
473.73 Hz: magnitude = 0.034186
495.26 Hz: magnitude = 0.005429
516.80 Hz: magnitude = 0.001862
538.33 Hz: magnitude = 0.000858
559.86 Hz: magnitude = 0.000466
581.40 Hz: magnitude = 0.000281
602.93 Hz: magnitude = 0.000182
624.46 Hz: magnitude = 0.000125

</pre>

We used a sin wave of 440Hz for the FFT analysis.
In the result, we can find that the bins that have the most highest magnitude are around 430Hz and 452Hz and which suggests the highest magnitude exists between these two bins(which is nearby 440Hz). 

The figure shows the spectrogram of the result.

![spectram_fft01.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/12cfccd9-d64c-4574-3304-c81291a59eec.jpeg)


##FFT()

This section instructs the user-function "FFT()"
The code of FFT() as follows,

<pre class="brush: c++;">
void FFT(float *wavedata, DSPSplitComplex *splitComplex)

{   
    int ffthalfsize = FFTSIZE/2;
    vDSP_Length log2n = log2(FFTSIZE);
    
    //setup fft
    FFTSetup fftsetup = vDSP_create_fftsetup(log2n, kFFTRadix2);
    
    //make window (fft size)
    float *window = malloc(sizeof(float)*FFTSIZE);
    //hanning window
    vDSP_hann_window(window, FFTSIZE, 0);
    //windowing
    vDSP_vmul(wavedata, 1, window, 1, wavedata, 1, FFTSIZE);

      
    //transform to complex
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );

    //FFT forward
    vDSP_fft_zrip(fftsetup, splitComplex, 1, log2n, FFT_FORWARD);

    //scaling
    float scale = 1.0/(FFTSIZE*2);
    vDSP_vsmul(splitComplex-&gt;realp, 1, &amp;scale, splitComplex-&gt;realp, 1, ffthalfsize);
    vDSP_vsmul(splitComplex-&gt;imagp, 1, &amp;scale, splitComplex-&gt;imagp, 1, ffthalfsize);

    //free
    vDSP_destroy_fftsetup(fftsetup);
    free(window);

}
</pre>

**FFT()** takes two argments, 1) float type of audio data buffer whose length is fft size, and **DSPSplitComplex型** data in which the results are stored as real and imaginary.


make Half FFT size

<pre class="brush: c++;">
int ffthalfsize = FFTSIZE/2;
vDSP_Length log2n = log2(FFTSIZE);
</pre>

Define setup of FFT for vDSP

<pre class="brush: c++;">
//setup fft
 FFTSetup fftsetup = vDSP_create_fftsetup(log2n, kFFTRadix2);
</pre>

The first argument is natural logarithm of FFT size, adn the second argument is the mode of FFT, in thise case I setup **kFFTRadix2 = 0**


Subsequently, generating a window functon. In thise case, I use hanning window.


<pre class="brush: c++;">
float *window = malloc(sizeof(float)*FFTSIZE);
vDSP_hann_window(window, FFTSIZE, 0); 
</pre>

the plot of a hanning window.

![spectram_fft02.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/8f1f9361-5336-79ef-62ad-35f80d4df986.jpeg)

multiply the generated hanning window with the audio data buffer.

<pre class="brush: c++;">
     vDSP_vmul(wavedata, 1, window, 1, wavedata, 1, FFTSIZE);
</pre>

The plot of the original audio data buffer

![spectram_fft03.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/8a733e99-84d0-eb87-247e-a4e7af1e91e8.jpeg)

The plot of the windowed audio data buffer


![spectram_fft04.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/f993eb15-3776-65da-a4ee-9456c8e8a4cf.jpeg)


Finally we can do FFT process!
Now we have to convert the windowed audio data buffer to **COMPLEX** data tyoes with which process is a bit tricky and very important to obtain the accurate result.

we use following function.

<pre class="brush: c++;">
    vDSP_ctoz( ( COMPLEX * )wavedata, 2, splitComplex, 1, ffthalfsize );
</pre>

We have to pay attention about the **Stride** arguments as they are set 2 for Complex wavedata and 1 for splitcomplex. It is because of that, audio data buffer (wavedata) has FFT size length and on the other hand, the length pof splitComplex is only half-fftsize. Thefore we need to set double number of splicomplex stride to the (Complex *)wavedata.

However, the process data amount should be ffthaldsize as it is set as 5th argument.

Now audio data buffer is converted to splitComplex. We can simply FFT splitComplex.

<pre class="brush: c++;">
     vDSP_fft_zrip(fftsetup, splitComplex, 1, log2n, FFT_FORWARD);
</pre>

Forward FFT is done by **FFT_FORWARD** as 5th argumtne and inverse FFT is done by **FFT_INVERSE**


We must not forget to scale the obtained splitComplex variables.
This process is important to get actual magnitude of the analyzed audio data and for the case you want to do inverse FFT afterwards.

<pre class="brush: c++;">
    float scale = 1.0/(FFTSIZE*2);
    vDSP_vsmul(splitComplex->realp, 1, &scale, splitComplex->realp, 1, ffthalfsize);
    vDSP_vsmul(splitComplex->imagp, 1, &scale, splitComplex->imagp, 1, ffthalfsize);
</pre>


Free all unneccesary data.
 

<pre class="brush: c++;">
    //free 
    vDSP_destroy_fftsetup(fftsetup);
    free(window);
</pre>


##FFT_free()						

also free splitComplex data

<pre class="brush: c++;">
void FFT_free(DSPSplitComplex *splitComplex)
{
    free(splitComplex->realp);
    free(splitComplex->imagp);
}
</pre>

##Square wave	


The Spectrogram of the square wave analyzed by above program.
		
![spectram_fft05.jpg](https://qiita-image-store.s3.amazonaws.com/0/46509/77336775-2682-9f5c-636b-04b9a23a8171.jpeg)