用c ++绘制频谱
Posted
技术标签:
【中文标题】用c ++绘制频谱【英文标题】:Plotting frequency spectrum with c++ 【发布时间】:2015-11-23 10:48:41 【问题描述】:请参阅此问题下方答案中的编辑内容。
我编写了一个脚本来用 c++ 绘制正弦信号的频谱。以下是步骤
-
应用汉宁窗
使用 fftw3 库应用 FFT
我有三张图:信号、信号何时乘以汉宁函数,以及频谱。频谱看起来不对。它应该在 50 Hz 处有一个峰值。任何建议将不胜感激。代码如下:
#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] 应该给我频谱。
以下是图表: 频谱和正弦时间窗口分别为:
【问题讨论】:
【参考方案1】:我的频率间隔完全错误。根据http://www.ni.com/white-paper/3995/en/#toc1; x 轴上的频率范围和分辨率取决于采样率和 N。频率轴上的最后一点应该是 Fs/2-Fs/N 和分辨率 dF=FS/N。所以我将脚本更改为:(因为频率分辨率为 Fs/N,当您增加采样数 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;
【讨论】:
记住t[i]
是没有意义的:你从不使用数组。此外,double t[N]
实际上不是有效的 C++。大小可变时使用std::vector<double>
,大小不变时使用const int N = 500
。
@MSalters ..当你说“大小是可变的”时,你能详细说明一下吗?向量的大小还是变量本身的大小?
我已经使用 t[i] 生成正弦波形
我的意思是当数组的大小不是运行时常数时,但两种解释实际上归结为相同。至于对波形使用t[i]
,您不使用那里的数组,而只是使用最后一个计算的元素。如果您将t[i]
替换为t[0]
,您的程序仍然有效!你可以在循环内声明double t = i * T
。
我一直在关注这个,但是使用了一个不受商业限制的库,我似乎没有在 out 数组上具有其他维度,使用仅标题的 FFT2D 库,它是真正的 dft 1 维函数似乎不会以相同的方式输出数据。任何想法如何将数组的就地更改映射到数据库?我复制了相同的正弦波函数,当我绘制它时,我得到看起来像噪声的值在 0 附近反弹。看不到波。【参考方案2】:
由于您的问题是关于 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;
您的代码可以使用更多 cmets,尤其是在循环或函数调用之前指定其输入值(目的)和/或返回值(结果)。
【讨论】:
【参考方案3】:我认为您可能没有足够的样本,特别是请参考此 Electronics.StackExhcange 帖子:https://electronics.stackexchange.com/q/12407/84272。
您要对 50 个样本进行采样,即 25 个 FFT 箱。您以 1000 Hz 的频率进行采样,因此每个 FFT 箱 1000 / 2 / 25 == 250 Hz。您的 bin 分辨率太低。
我认为您需要降低采样频率或增加采样数。
【讨论】:
以上是关于用c ++绘制频谱的主要内容,如果未能解决你的问题,请参考以下文章
Merry Christmas!一起用C语言绘制一个动态的圣诞树吧
Merry Christmas!一起用C语言绘制一个动态的圣诞树吧