Jack
Jack

Reputation: 753

Spectrogram with C++

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

  1. I divided the real sinusoidal signal into B blocks
  2. Applied Hanning window to each block (I assumed there is no overlap). This should give me the inputs for the fft.
  3. 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
  4. Plotted time tt[k] vs v[k][i]

This gives me the following plot:

enter image description here

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;
}

Upvotes: 0

Views: 1686

Answers (1)

Jack
Jack

Reputation: 753

enter image description here

#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;
}

Upvotes: 3

Related Questions