Please see the Edits in the answer below this question.
I have written a script to plot the frequency spectrum of a sinusoidal signal with c++. Here are the steps
- Applying Hanning window
- Apply FFT using fftw3 library
I have three graphs: Signal, Signal when is multiplied to Hanning function, and the frequency spectrum. The frequency spectrum looks wrong. It should have a peak at 50 Hz. Any suggestion would be appreciated. Here is the 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;
double y;
int N=50;
double Fs=1000;//sampling frequency
double T=1/Fs;//sample time
double f=50;//frequency
double *in;
fftw_complex *out;
double t[N];//time vector
double ff[N];
fftw_plan plan_forward;
in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
for (int i=0; i< N;i++)
{
t[i]=i*T;
ff[i]=1/t[i];
in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
in[i] = multiplier * in[i];
}
plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );
fftw_execute ( plan_forward );
double v[N];
for (int i = 0; i < N; i++)
{
v[i]=20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])/N/2);//Here I have calculated the y axis of the spectrum in dB
}
fstream myfile;
myfile.open("example2.txt",fstream::out);
myfile << "plot '-' using 1:2" << std::endl;
for(i = 0; i < N; ++i)
{
myfile << ff[i]<< " " << v[i]<< std::endl;
}
myfile.close();
fftw_destroy_plan ( plan_forward );
fftw_free ( in );
fftw_free ( out );
return 0;
}
I have to add that I have plotted the graphs using gnuplot after inserting the results into example2.txt. So ff[i] vs v[i] should give me the frequency spectrum.
Here are the plots: Frequency Spectrum and Sinusoidal time Window respectively:
I think you may not have enough samples, particularly, reference this Electronics.StackExhcange post: https://electronics.stackexchange.com/q/12407/84272.
You're sampling for 50 samples, so 25 FFT bins. You're sampling at 1000 Hz, so 1000 / 2 / 25 == 250 Hz per FFT bins. Your bin resolution is too low.
I think you need to lower the sampling frequency or increase the number of samples.
Since your question in on SO, your code could use some indentation and style improvement to make it easier to read.
Your code can use more comments, especially before loops or function calls to specify their input value (purpose) and/or returning value (result).