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.
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