FFT of wav file using C++ / Python

Question

I have been googling extensively and could plot the FFT of my wav file using Python but am unable to do so for C++, which I originally had to do. I downloaded and linked the FFTW and LIBSND to Visual C++. Though I am not understanding which functions to use in the library and what to send to compute the same kind of result as I did in Python.

Question : I basically first am not understanding, how to get the arrays of Amplitude, Frequency and then I'll go ahead and plot them.

My C++ code for the same is : `

//.........................np.fft.fft as in python
p = fftw_plan_dft_1d(num_items, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); 
fftw_execute(p);

//According to Nyquist its 1/2
for (int i = 0; i < num_items / 2; ++i) {
            printf("real=%f ",out[i][0]);
            printf("img=%f ",out[i][1]);
        }

//Amplitude
float *amp;
amp = (float *)malloc(sizeof(float)*(num_items / 2));
for (int i = 0; i < num_items/2; ++i) {
    amp[i] = sqrt( pow(out[i][0],2) + pow(out[i][1], 2));
}`

Here's the python code for it.

import sys
import numpy as np
from scipy.io.wavfile import read
from matplotlib import pyplot as plt


def do_fft(received_wave, Fs=44100):
    """
    :param received_wave: wave file data.
    :param Fs: Sampling Rate, default = 44100
    :return: [Frequency, Amplitude]
    """

    # Calculating the fft coeff and amp sqrt(x^2+y^2)
    fft_coeff   = np.fft.fft(received_wave)
    Amp         = np.sqrt(np.abs(fft_coeff))


    print "FFT_coeff: ",fft_coeff
    print "Amp: ",Amp 

    # calulating size of recieved wave data and creating a freq array based on sampling freq Fs and size
    size1=len(received_wave)
    freq=np.linspace(0,Fs,size1)

    print "Length of recieved wave: ",size1;

    # Taking only half sample based on Nyquist-Shannon sampling theorem for ampiltude and frequency
    # https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem
    Amplitude  = Amp[0:int(size1/2)]
    Frequency = freq[0:int(size1/2)]


    print "\nAmplitude : ", Amplitude
    print "\nFreq : ", Frequency
    # This shorts the  index of the array in acending order
    idx = np.argsort(Amplitude)
    # freq1 is the maximum freq freq2 second maximum and so on
    freq1 = ((idx[-1]) / float(size1)) * Fs  
    freq2 = ((idx[-2]) / float(size1)) * Fs
    freq3 = ((idx[-3]) / float(size1)) * Fs

    return Amplitude, Frequency, freq1, freq2, freq3


def read_from_file(file_location):
    """
    Read file ad return audio data
    :param file_location: location of file.
    :return: audio data
    """

    data = read(file_location)


    # as scipy read function return two array [sample_rate_of_file, [audio_chunks]]
    sample_rate, audio_data = data
    print "Data: " ,data
    for i in range(len(audio_data)):
     #   print audio_data[i]
         pass

    print i
    return sample_rate, audio_data

def plot_fft(audio_file):

    # read audio chunks from audio file
    sample_rate, audio_data = read_from_file(audio_file)

    # call do_fft() function to get fft ( frequency and amplitude)
    Amplitude, Frequency, freq1, freq2, freq3 = do_fft(received_wave=audio_data, Fs=sample_rate)

    # plot fft
    plt.title("FFT heigest : {}, second_heigest : {}".format(freq1, freq2))
    plt.plot(Frequency, Amplitude)
    plt.show()
    plt.close()

    return True


if __name__ == '__main__':
    file = "hellotrill.wav"
    plot_fft(file)

I'm not understanding as the arrays I'm getting in Python and C++ completely differ.

score 3 · Answer 1 · answered Apr 15 '17 at 11:35

3

If you want to plot something you need to use a library that can draw graphics or generate data for another program that can.

There are many options - a few you could use would be:

The C++ standard library does not provide any graphics capabilities out-of-the-box.

answered Apr 15 '17 at 11:35

Jesper Juhl

30,449
3
47
70

Thanks for the answer. Though I was looking more along the lines of , how to compute FFT and get the amplitude and frequency, freq1,freq2,fre3 as shown in python. I will also edit the question again. – Shivansh Jagga Apr 15 '17 at 18:08

Shivansh Jagga · Accepted Answer · 2017-04-19T15:24:56.927

I converted the code and successfully plotted it using Gnuplot. The code is given below and was done in Visual Studio 2017 (hence the stdafx.h header file)

#include "stdafx.h"
#include<fftw3.h>
#include<sndfile.h>

#include<fstream>
#include<vector>
#include<math.h>
#include <algorithm>
#include<iostream>
using namespace std;

#define file_path "F:/Shivansh Work/University work/VIT/Assignment2/hellotrill.wav"


//Read http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs

void plot_fft(float *amp, float *freq, float *freq2,int num_items, int Fs) {

fstream amp_freq;
amp_freq.open("fft_plot.txt", fstream::out);
for (int i = 0; i < num_items/2; ++i)
{
    amp_freq << freq2[i] << " " << amp[i] << std::endl;
}
amp_freq.close();


float *idx; // amplitude array for sorting
idx= (float *)malloc(sizeof(float)*(num_items / 2));
for (int i = 0; i < num_items / 2; ++i) {
    idx[i] = amp[i] ;   
}

int size= num_items / 2;
sort(idx, idx + size);
cout << idx[size - 1] << " " << amp[size - 1];
//NOTE: np.argsort returns the indices of sorted array, but not the values itself
float fre1, fre2, fre3;
fre1 = Fs*idx[size - 1]/ (float)num_items;
fre2 = Fs*idx[size - 2]/ (float)num_items;
fre3 = Fs*idx[size - 3]/ (float)num_items;

printf("\n\nHighest frequencies: %.5f, %.5f, %.5f", fre1, fre2, fre3);
printf("\n[NOTE: In given python code: np.argsort returns the indices of sorted array, but not the values itself]");
printf("\n\nThe amplitude and frequency have been written in the file fft_plot.txt");

}

int main() {

char        *infilename;
SNDFILE     *file = NULL;

SF_INFO     sfinfo;
int num_channels;
int num, num_items;
double *buf;
int frame, samplerate, ch;
int i, j;

FILE        *outfile = NULL;

//Read the file, into buffer
file = sf_open(file_path, SFM_READ, &sfinfo);


/* Print some of the info, and figure out how much data to read. */
frame = sfinfo.frames;
samplerate = sfinfo.samplerate;
ch = sfinfo.channels;
printf("frames=%d\n", frame);
printf("samplerate=%d\n", samplerate);
printf("channels=%d\n", ch);
num_items = frame * ch;
printf("num_items=%d\n", num_items);



//Allocate space for the data to be read, then read it
buf = (double *)malloc(num_items * sizeof(double));
num = sf_read_double(file, buf, num_items);

sf_close(file);
printf("Read %d items\n", num);


/*initialize FFT parameters*/
fftw_complex *in, *out;
fftw_plan p;
/*Do fft to wav data*/
in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * num_items);
out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * num_items);
for (int i = 0; i < num_items; i++) {
    in[i][0] = buf[i];
    in[i][1] = 0;
}
//.........................np.fft.fft as in python (OUT stores fft_coeff)
p = fftw_plan_dft_1d(num_items, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);     //1D Complex DFT, FFTW_FORWARD & BACKWARD just give direction and have particular values
fftw_execute(p);


/*  //According to Nyquist its 1/2
for (int i = 0; i < num_items / 2; ++i) {
            printf("%f+",out[i][0]);
            printf("%fj  ",out[i][1]);
        }
*/

//...............................Amplitude
float *amp;
amp = (float *)malloc(sizeof(float)*(num_items / 2));
for (int i = 0; i < num_items/2; ++i) {
    amp[i] = sqrt (sqrt( pow(out[i][0],2) + pow(out[i][1], 2)));        //2 sqrt since np.sqrt( np.abs() )
}
fftw_destroy_plan(p);


//...............................Frequency
float *freq;
freq = (float *)malloc(sizeof(float)*(num_items/2));
int size = samplerate / num_items;

double *samples;
samples = (double *)malloc(sizeof(double)*samplerate);     //Multiplying by sample rate cuz of np.linspace, goes from 0-samplerate
sf_read_double(file, samples, samplerate);

fftw_complex* out2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * num_items);
fftw_plan plan;
plan = fftw_plan_dft_r2c_1d(num_items, samples, out2, FFTW_ESTIMATE);   //out2 imaginary parts are all 0, can read http://stackoverflow.com/questions/4364823/how-do-i-obtain-the-frequencies-of-each-value-in-an-fft
fftw_execute(plan);

for (int i = 0; i<num_items/2; i++) 
{
    freq[i] = sqrt(pow(out2[i*size][0],2) + pow(out2[i*size][1],2));//np.linspace (0-Fs ,in size1 increments), also can read http://stackoverflow.com/questions/4364823/how-do-i-obtain-the-frequencies-of-each-value-in-an-fft
}
//NOTE: In np.linspace(0,44100,29757) -> a normal array is created with numbers. Not actual frequency. 
//But here actual frequency is being created

fftw_destroy_plan(plan);
fftw_free(out2); fftw_free(in); fftw_free(out);


//................................Function for frequency according to python program
float *freq2;
freq2 = (float *)malloc(sizeof(float)*(num_items / 2));

float size2 = ( float)samplerate / (float)num_items;
for (int i = 0; i<= num_items / 2; i++)
{
    freq2[i] = i*size2;
}

//..................................Plotting
plot_fft(amp, freq, freq2, num_items,samplerate);

return 0;

}

The array was turned into a text file as shown in the plot_fft function. And then was plotted in gnuplot easily.

PLOTTING Function set style line 1 lc rgb '#0060ad' lt 1 lw 2 pt 7 ps 1.5 #--- blue plot 'fft_plot.dat' with linespoints ls 1

FFT of wav file using C++ / Python

2 Answers2

Linked