spectrum and oscillation frequency of a vector

[jagc944]

I'm trying to get the spectrum and oscillation frequency of a pixel vector (motion of a point in a video).

The vector data (file.txt) is organized as follows:

21.000000
24.000000
25.000000
25.000000
21.000000
22.000000
24.000000
25.000000
24.000000
22.000000
114.000000
75.000000
60.000000
52.000000
61.000000
66.000000
78.000000
86.000000
74.000000
59.000000
47.000000
58.000000
60.000000
81.000000
85.000000
81.000000
70.000000
58.000000
59.000000
56.000000
61.000000
77.000000
88.000000
82.000000
79.000000
75.000000
75.000000
75.000000
75.000000
76.000000
82.000000
82.000000
85.000000
85.000000
78.000000
74.000000
77.000000
80.000000
80.000000
81.000000

I have read a bit about the possibility of doing this with the fftw3 library, but I am new to handling c ++ and this library.

I need help with codes for obtaining the spectrum (db) and the frequency (Hz) of my vector.

Thank you very much for your help

[lastchance]

Not sure exactly what you want - there isn't a single "frequency" in your system: it's a combination of frequencies which make up a "spectrum".

The code below will give you some (discrete) Fourier transforms. The columns in the output file are:
- frequency (as a multiple of the base frequency: 1/sample length)
- real part of discrete Fourier transform
- imaginary part of discrete Fourier transform
- square of the magnitude of the Fourier transform ("spectral energy")

The sample points you have given are overwhelmed by the low-frequency component, but I assume it's not the whole sample.

#include <iostream>
#include <fstream>
#include <iomanip>
#include <complex>
#include <cmath>
using namespace std;

const double pi = 3.14159265358979323846;

void FFT ( int N, complex<double> *f     , complex<double> *ftilde );
void iFFT( int N, complex<double> *ftilde, complex<double> *f      );
void DFT ( int N, complex<double> *f     , complex<double> *ftilde );
void iDFT( int N, complex<double> *ftilde, complex<double> *f      );
void writeFile( string filename, int N, complex<double> *z );


//======================================================================


int main()
{
   string infile = "file.txt";
   string oufile = "fourier.out";
   string line;
   double re, im;

   // First pass - count data (assumes one per line)
   int N = 0;
   ifstream in( infile.c_str() );
   while( getline( in, line ) ) N++;
   cout << "Number of points read from file: " << N << endl;
  
   // Reposition at start
   in.close();
   in.open( infile.c_str() );

   // Reread, filling f
   complex<double> *f = new complex<double>[N];
   for ( int i = 0; i < N; i++ )
   {
      in >> re;   im = 0.0;        // adjust if you have complex data
      f[i] = complex<double>( re, im );
   }
   in.close();

   // Fourier transform
   complex<double> *ftilde = new complex<double>[N];
   FFT( N, f, ftilde );

   // Output to file
   writeFile( oufile, N, ftilde );

   // Tidy up
   delete [] f;
   delete [] ftilde;
}


//======================================================================


void FFT( int N, complex<double> *f, complex<double> *ftilde )       // fast fourier transform
{
   if ( N % 2 == 0 )                                                 // recursion if a multiple of 2
   { 
      int N2 = N / 2;
      complex<double> *even = new complex<double>[N2];
      complex<double> *odd  = new complex<double>[N2];
      complex<double> *g    = new complex<double>[N2];
      complex<double> *h    = new complex<double>[N2];

      for ( int m = 0; m < N2; m++ )
      {
         even[m] = f[2*m  ];
         odd [m] = f[2*m+1];
      }

      FFT( N2, even, g );
      FFT( N2, odd , h );

      complex<double> w = polar( 1.0, -2.0 * pi / N );
      complex<double> wn = complex<double>( 1.0 );
      for ( int n = 0; n < N2; n++ )
      {
         int nplus = n + N2;
         ftilde[n]     = g[n] + wn * h[n];
         ftilde[nplus] = g[n] - wn * h[n];
         wn *= w;
      }

      delete [] even;
      delete [] odd ;
      delete [] g   ;
      delete [] h   ;
   }

   else                                                              // otherwise just do a DFT
   {
      DFT( N, f, ftilde );
   }
}


//======================================================================


void iFFT( int N, complex<double> *ftilde, complex<double> *f )      // inverse fast fourier transform
{
   complex<double> * ftildeConjugate = new complex<double>[N];

   for ( int m = 0; m < N; m++ ) ftildeConjugate[m] = conj( ftilde[m] );

   FFT( N, ftildeConjugate, f );

   for ( int m = 0; m < N; m++ ) f[m] = conj( f[m] ) / (double) N;

   delete [] ftildeConjugate;
}


//======================================================================


void DFT( int N, complex<double> *f, complex<double> *ftilde )       // discrete Fourier transform
{
   for ( int n = 0; n < N; n++ )
   {
      ftilde[n] = complex<double>( 0.0 );
      complex<double> w = polar( 1.0, -2.0 * pi * n / N );
      complex<double> wm( 1.0 );
      for ( int m = 0; m < N; m++ )
      {
         ftilde[n] += f[m] * wm;
         wm *= w;
      }
   }
}


//======================================================================


void iDFT( int N, complex<double> *ftilde, complex<double> *f )      // inverse discrete Fourier transform
{
   complex<double> * ftildeConjugate = new complex<double>[N];

   for ( int m = 0; m < N; m++ ) ftildeConjugate[m] = conj( ftilde[m] );
   DFT( N, ftildeConjugate, f );
   for ( int m = 0; m < N; m++ ) f[m] = conj( f[m] ) / (double) N;

   delete [] ftildeConjugate;
}


//======================================================================


void writeFile( string filename, int N, complex<double> *z )
{
   double freq, re, im, magsq;

   ofstream out( filename.c_str() );

   for ( int m = 0; m < N; m++ ) 
   {
      re = real( z[m] );
      im = imag( z[m] );
      magsq = re * re + im * im;
      out << right
          << setw(7) << m << "  "
          << scientific
          << setprecision(6) << setw(15) << re << "  "
          << setprecision(6) << setw(15) << im << "  "
          << setprecision(6) << setw(15) << magsq << endl;
   }

   out.close();
}


//======================================================================