Spectrogram with C++

Question

I am trying to generate a spectrogram of a sinusoidal signal with C++. To generate the spectrogram:

I divided the real sinusoidal signal into B blocks
Applied Hanning window to each block (I assumed there is no overlap). This should give me the inputs for the fft.
Applied fft to each block and calculated the magnitude from the frequency coefficients u[i][0] and u[i][1] and put it into v[k][i] where k is the number of blocks and i the samples of each window
Plotted time tt[k] vs v[k][i]

This gives me the following plot:

However, that plot does not look right.

Any suggestions to make it work? Here is my code:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
int main()
{
  int i;
  int N=500;//Number of points acquired inside the window
  double Fs=200;//sampling frequency
  int windowsize=100;//Number of samples in each block
  double  T=1/Fs;//sample time 
  double f=50;//frequency(Hz)
  double *in;
  fftw_complex *out;
  double t[N];//time vector 
  double tt[5];
  double  v [5][249];
  fftw_plan plan_forward;
  in = (double*) fftw_malloc(sizeof(double) * N);
  out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));

  for (int i=0; i<= N;i++)
  {
    t[i]=i*T;

    in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
  }

  plan_forward = fftw_plan_dft_r2c_1d (windowsize, in, out, FFTW_ESTIMATE );

  for (int k=0; k< 5;k++){ 
    for (int i = 0; i<windowsize; i++){
      double multiplier = 0.5 * (1 - cos(2*M_PI*i/(windowsize-1)));//Hanning  Window
      in[i] = multiplier * in[i+k*windowsize];
      fftw_execute ( plan_forward );
      v[k][i]=(20*log10(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])))/N;//The magnitude in dB
    }
  }

  for (int k=0; k< 5;k++){//Center time for each block

    tt[k]=(2*k+1)*T*(windowsize/2);

  }

  fstream myfile;
  myfile.open("example2.txt",fstream::out);

  myfile << "plot '-' using 1:2" << std::endl;

  for (int k=0; k< 5;k++){ 
    myfile << v[k][i]<< " " << tt[k]<< std::endl;
  }
  myfile.close();
  fftw_destroy_plan ( plan_forward );
  fftw_free ( in );
  fftw_free ( out );
  return 0;
}

I'll try to post a more complete answer later, but quickly the most obvious issue is that you are only writing tt (instead of both the time & the data) to file, then try to plot it as a x vs. y function (whereas you really have 2D data as a function of time, so you'd need a 3D plot). Also, fftw_execute should be called on each block of data, not for each individual sample fed into your input buffer. — SleuthEye
That was a typo. I tried to see the results of tt and forgot to put it back to tt vs v but the plot is for myfile << v[k][i]<< " " << tt[k]<< std::endl; — Jack
typo aside, 2D plot is still not the right plot type to show a spectrogram where you have intensity in dB for all times and frequency bins (ie. 3D). The intensity is often represented by a color mapping so it can be seen on a 2D plane. — SleuthEye

Jack Jack · Accepted Answer · 2015-09-03T19:10:11

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int main()
{
int i;
int N=500;
double Fs=200;//sampling frequency
int windowsize=100;//Number of points acquired inside the window
double dF=Fs/N;
double  T=1/Fs;//sample time 
double f=50;//frequency
double *in;
fftw_complex *out;
double t[N];//time vector 
double ff[N];
double tt[5];
double  v [5][249];
fftw_plan plan_forward;
in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));

for (int i=0; i<= N;i++)
{ t[i]=i*T;
in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
 }

 plan_forward = fftw_plan_dft_r2c_1d (windowsize, in, out, FFTW_ESTIMATE );

for (int k=0; k< 5;k++){ 
  for (int i = 0; i<windowsize; i++){
  double multiplier = 0.5 * (1 - cos(2*M_PI*i/(windowsize-1)));//Hanning Window
 in[i] = multiplier * in[i+k*windowsize];
                                     }
fftw_execute(plan_forward);
for (int j = 0; j <= ((windowsize/2)-1); j++){
v[k][j] = 2*(out[j][0] * out[j][0] + out[j][1] * out[j][1])/(N*N);
v[k][0]   *= 0.5;
 v[k][N/2] *= 0.5;

for (int j = 0; j <= windowsize/2; j++){
   v[k][j] = 10 * log10(v[k][j] + 1e-5);
                                        }
                  }

for (int j = 0; j <= ((windowsize/2)-1); j++)
    {ff[j]=Fs*j/windowsize;
     }

//Center time for each block

for (int k=0; k< 5;k++){

tt[k]=(2*k+1)*T*(windowsize/2);
                        }

double b[6];
fstream myfile;

myfile.open("data.txt",fstream::out);

myfile << "plot '-' using 1:2" << std::endl;

for (int k=0; k< 6;k++) {
                    b[0]=5;
                    b[k+1]=tt[k];
                     myfile <<b[k] << "\t";
                    }
 myfile<< std::endl;

 for (int j = 0; j <= windowsize/2; j++){  myfile << ff[j]<< "\t";

   for (int k=0; k< 5;k++){ myfile << v[k][j]<< "\t";

                          }
 myfile<< std::endl;
                                    }
 myfile.close();

 fftw_destroy_plan ( plan_forward );
 fftw_free ( in );
 fftw_free ( out );
 return 0;
}

Spectrogram with C++

1 Answers