Построение частотного спектра с переполнением стека
Пожалуйста, смотрите Редактирование в ответе ниже этого вопроса.
Я написал скрипт для построения частотного спектра синусоидального сигнала с помощью c ++. Вот шаги
- Применение окна Хеннинга
- Применить БПФ с помощью библиотеки fftw3
У меня есть три графика: сигнал, сигнал при умножении на функцию Хеннинга и частотный спектр. Частотный спектр выглядит неправильно. Он должен иметь пик при 50 Гц. Любое предложение будет оценено. Вот код:
#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;
}
Я должен добавить, что я построил графики с помощью gnuplot после вставки результатов в example2.txt. Поэтому ff [i] vs v [i] должен дать мне частотный спектр.
Вот графики: 
Частотный спектр и синусоидальное временное окно соответственно:
Решение
Мои частотные интервалы были совершенно неверными. В соответствии с http://www.ni.com/white-paper/3995/en/#toc1; частотный диапазон и разрешение на Икс-ось зависит от частоты дискретизации и N. Последняя точка на оси частот должна быть Fs / 2-Fs / N и разрешение DF = FS / N.Поэтому я изменил свой сценарий на: (так как разрешение по частоте Fs / N по мере увеличения количества smaples N (или уменьшить частоту дискретизации Fs) вы получите меньшее разрешение по частоте и лучшие результаты.)
#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=550;//Number of points acquired inside the window
double Fs=200;//sampling frequency
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];
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;
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];
}
for (int i=0; i<= ((N/2)-1);i++)
{ff[i]=Fs*i/N;
}
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/2)-1); i++)
{
v[i]=(20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])))/N; //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/2)-1); i++)
{
myfile << ff[i]<< " " << v[i]<< std::endl;
}
myfile.close();
fftw_destroy_plan ( plan_forward );
fftw_free ( in );
fftw_free ( out );
return 0;
}
Другие решения
Я думаю, что вам может не хватить образцов, в частности, сослаться на этот пост Electronics.StackExhcange: https://electronics.stackexchange.com/q/12407/84272.
Вы отбираете 50 образцов, так что 25 бин FFT. Вы производите выборку с частотой 1000 Гц, поэтому 1000/2/25 == 250 Гц на бункеры FFT. Разрешение вашего бина слишком низкое.
Я думаю, вам нужно уменьшить частоту дискретизации или увеличить количество выборок.
Поскольку ваш вопрос касается SO, ваш код может использовать некоторые отступы и улучшения стиля, чтобы его было легче читать.
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
int main(){
// use meaningful names for all the variables
int i;
double y;
int N = 550; // number of points acquired inside the window
double Fs = 200; // sampling frequency
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];
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;
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];
}
for(int i = 0; i <= ((N/2)-1); i++){
ff[i] = (Fs * i) / N;
}
plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
fftw_execute(plan_forward);
double v[N];
// Here I have calculated the y axis of the spectrum in dB
for(int i = 0; i <= ((N/2)-1); i++){
v[i] = (20 * log(sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]))) / N;
}
fstream myfile;
myfile.open("example2.txt", fstream::out);
myfile << "plot '-' using 1:2" << std::endl;
for(i = 0; i < ((N/2)-1); i++){
myfile << ff[i] << " " << v[i] << std::endl;
}
myfile.close();
fftw_destroy_plan(plan_forward);
fftw_free(in);
fftw_free(out);
return 0;
}
Ваш код может использовать больше комментариев, особенно перед циклами или вызовами функций, чтобы указать их входное значение (цель) и / или возвращаемое значение (результат).