Search code examples
c++fftw

calculate the precise, non integer frequencies from time series in (FFTW)

I want to calculate the frequency of time series precisely with at least 3 decimal value. This is a simple example that calculates the frequency of integer values.

#include <fftw3.h>
#include <cstdio>
#include <cmath>
#include <iostream>
#include <fstream>

#define REAL 0
#define IMAG 1
#define NUM_POINTS 1024


void acquire_signal(double *signal, double *theta) {
    /* Generate two sine waves of different frequencies and
     * amplitudes.
     */

    int i;
    for (i = 0; i < NUM_POINTS; ++i) {
        theta[i] = (double)i / (double)NUM_POINTS;
        signal[i] = 1.0*sin(50.0 * 2.0 * M_PI * theta[i]) +
                    0.5*sin(80.0 * 2.0 * M_PI * theta[i]);
    }
}

int main() {
    unsigned flags{0};
    double *theta  = new double[NUM_POINTS];
    double *signal = new double[NUM_POINTS];

    fftw_complex result[NUM_POINTS/2+1];

    fftw_plan plan = fftw_plan_dft_r2c_1d(NUM_POINTS,
                                         signal,
                                         result,
                                         flags);
    acquire_signal(signal,theta);
    fftw_execute(plan);

    //save signal and result
    std::ofstream f1,f2;
    f1.open ("signal.txt");
    for (int i=0; i<NUM_POINTS; i++){
        f1 <<theta[i]<<" "<<signal[i]<<"\n";
    }

    f1.close();
    f2.open("result.txt");

    for (int i=0; i<NUM_POINTS/2; i++){
        double yf = 2.0/(double)(NUM_POINTS)* sqrt(result[i][REAL]*result[i][REAL]+ result[i][IMAG]* result[i][IMAG]);
        f2<< (double)i << " "<<yf <<"\n";
    }
    f2.close();
    fftw_destroy_plan(plan);
    delete[] signal,theta;

    return 0;
}

enter image description here But how should I change the code if I have 

signal = 1.0*sin(50.350 * 2.0 * M_PI * theta[i]) +
         0.5*sin(80.455 * 2.0 * M_PI * theta[i]);

Is it appropriate to change the units of time and frequency? for example time in 1000*s and frequency in kHz?

Solution

  • Just changing the numbers in sin will shift your lines from 50 and 80 to 50.350 and 80.455 Hz, and assuming you have 1024 lines by 1024 Hz. But you still have 1Hz resolution. You need more lines (x1000) by the same sampling frequency to get bigger resolution.

    For example if you want 1/4 Hz resolution you need 4x more samples so by 1024 Hz sample rate you need fs * 4 samples:

    ...

#define NUM_POINTS (1024 * 4)
double fs = 1024;  // Sample rate in Hz

void acquire_signal(double *signal, double *theta) {
    /* Generate two sine waves of different frequencies and
     * amplitudes.
     */

    int i;
    for (i = 0; i < NUM_POINTS; ++i) {
        theta[i] = (double)i / (double)fs;
        signal[i] = 1.0*sin(50.0 * 2.0 * M_PI * theta[i]) +
                    0.5*sin(80.0 * 2.0 * M_PI * theta[i]);
    }
}
....
for (int i=0; i< (NUM_POINTS/2 + 1) ; i++){
    double yf = 2.0/(double)(NUM_POINTS)* sqrt(result[i][REAL]*result[i][REAL]+ result[i][IMAG]* result[i][IMAG]);
    f2 << (double)i * fs / ( NUM_POINTS ) << " "<<yf <<"\n";
}

    0 2.90715e-16

    0.25 1.19539e-16

    0.5 2.15565e-16

    0.75 2.88629e-16

    1 3.05084e-16

    1.25 3.864e-16

    ...

    49.75 9.47968e-16

    50 1

    50.25 1.12861e-15

    50.5 4.95946e-16

    50.75 6.9016e-16

    ...