使用C ++绘制频谱 [英] Plotting frequency spectrum with c++
问题描述
请在此问题下方的答案中查看编辑".
Please see the Edits in the answer below this question.
我写了一个脚本，用c ++绘制正弦信号的频谱.这是步骤
I have written a script to plot the frequency spectrum of a sinusoidal signal with c++. Here are the steps
- 应用汉宁窗口
- 使用fftw3库应用FFT
我有三个图表:信号，信号何时乘以汉宁函数以及频谱.频谱看起来不对.它应该在50 Hz处有一个峰值.任何建议，将不胜感激.这是代码:
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;
}
我必须补充一点，在将结果插入example2.txt之后，我已经使用gnuplot绘制了图表. ff [i] vs v [i]应该给我频谱.
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:
推荐答案
我的频率间隔完全错误.根据 http://www.ni.com/white-paper/3995/en/#toc1 ; x 轴上的频率范围和分辨率取决于采样率和 N .频率轴上的最后一个点应该是 Fs/2-Fs/N 和分辨率 dF = FS/N .因此，我将脚本更改为:(由于频率分辨率为 Fs/N ，这是因为您增加了样本的数量 N (或降低采样频率 Fs )，从而获得了较小的频率分辨率和更好的结果. )
My Frequency intervals were completely wrong. According to http://www.ni.com/white-paper/3995/en/#toc1; the frequency range and resolution on the x-axis depend on sampling rate and N. The last point on the frequency axis should be Fs/2-Fs/N and the resolution dF=FS/N.So I have changed my script to: (since frequency resolution is Fs/N as you increase the number of smaples N (or decrease sampling frequency Fs) you get smaller frequency resolution and better results.)
#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;
}
这篇关于使用C ++绘制频谱的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
