有时候需要反复验证一点点数字信号处理的程序,因为简单,通常都是单个文件进行编译。如果只执行一次,将数据保存到文件里,然后使用电子表格打开,可以绘制曲线。但是如果需要反复的验证调整编写函数,那么能简单的绘制曲线就很有必要了。
因为是出于验证目的,速度就不是很必要,能验证结果就好。
据说gnuplot是能媲美Matlab的绘图软件。能不能媲美不重要,重要的是它小巧,但功能完备,简单易用。
- 先提供一个简易的gnuplot绘制曲线的封装接口。
// File Name: MySimplePlot.hpp
#ifndef _MY_SIMPLE_PLOT_H_
#define _MY_SIMPLE_PLOT_H_
#ifdef _WIN32
#include <windows.h>
#define popen _popen
#define pclose _pclose
#else
#include <unistd.h>
#endif
#include <vector>
#include <string>
#include <stdexcept>
#include <stdio.h>
// gnuplot的GUI有wxWidgets和Qt两种
#define USE_WX
//#define USE_QT
class XYChart
{
public:
XYChart();
~XYChart();
// 清除所有数据
void clear();
// 设置画图窗口的标题和XY坐标的标签
void title(const char* title, const char* xlabel=NULL, const char* ylabel=NULL );
// 添加数据
void add(const char* title,int length, double* y, double* x = NULL);
// 绘制直线
void plot();
private:
XYChart(const XYChart&);
XYChart& operator=(const XYChart& src);
void * operator new(size_t size);
FILE *gnu_pipe;
bool init_ok;
std::vector<std::vector<double>> Xs;
std::vector<std::vector<double>> Ys;
std::vector<std::string> Ts;
std::vector<std::string> Fs;
};
XYChart::XYChart()
{
gnu_pipe = popen("gnuplot", "w");
//gnu_pipe = popen("\"C:\\Program Files\\gnuplot\\bin\\gnuplot.exe\"", "w");
if(gnu_pipe)
{
init_ok = true;
#ifdef USE_QT
fprintf(gnu_pipe, "set term qt 0 title \"XYPlot\"\n");
#endif
#ifdef USE_WX
fprintf(gnu_pipe, "set term wx 0 title \"XYPlot\"\n");
#endif
}
else throw std::runtime_error("XYChart() cannot open gnuplot");
}
XYChart::~XYChart()
{
if(init_ok)
{
fprintf(gnu_pipe, "exit\n");
pclose(gnu_pipe);
clear();
}
}
void XYChart::title( const char* title, const char* xlabel, const char* ylabel )
{
if(init_ok)
{
fprintf(gnu_pipe, "set title \"%s\"\n", title);
if(xlabel) fprintf(gnu_pipe, "set xlabel \"%s\"\n", xlabel);
if(ylabel) fprintf(gnu_pipe, "set ylabel \"%s\"\n", ylabel);
}
}
void XYChart::add(const char* title,int N, double* y, double* x )
{
// 添加数据,就是简单的将数据写入文件
if(!title)
throw std::invalid_argument("XYChart::add() the title is NULL");
if(N <= 0)
throw std::invalid_argument("XYChart::add() the data length N <= 0");
std::vector<double> X(N), Y(N);
for(int i = 0; i < N; i++) Y[i] = y[i];
if(x != NULL)
{
for(int i = 0; i < N; i++) X[i] = x[i];
}
else
{
for(int i = 0; i < N; i++) X[i] = i;
}
std::string T(title);
#ifdef _WIN32
std::string filename = "D:\\\\";
#else
std::string filename = "./";
#endif
filename += T;
filename += ".lqbz";
Fs.push_back(filename);
Xs.push_back(X);
Ys.push_back(Y);
Ts.push_back(T);
FILE* pf = fopen(filename.c_str(), "w");
if(pf)
{
for(int n = 0; n < N; n++)
fprintf(pf, "%lf %lf \n", X[n], Y[n]);
fclose(pf);
}
else throw std::runtime_error("XYChart::add() cannot open file for writing data");
}
void XYChart::clear()
{
Xs.clear();Ys.clear();Ts.clear();
int count = Fs.size();
for(int i = 0; i < count; i++) remove(Fs[i].c_str());
Fs.clear();
}
void XYChart::plot()
{
// 直接通过管道调用gnuplot的指令画图
// 可以很方便的修改成其他需要的命令
int array_count = Xs.size();
if(array_count == 0) return;
if(init_ok)
{
fprintf(gnu_pipe, "clear\n");
fprintf(gnu_pipe, "unset label\n");
fprintf(gnu_pipe, "set autoscale xy\n");
fprintf(gnu_pipe, "plot \"%s\" using 1:2 title \"%s\" with lp lw 2 pt 0\n", Fs[0].c_str(), Ts[0].c_str());
for(int i = 1; i < array_count; i++)
{
fprintf(gnu_pipe, "replot \"%s\" using 1:2 title \"%s\" with lp lw 2 pt 0\n", Fs[i].c_str(), Ts[i].c_str());
}
fprintf(gnu_pipe, "\n\n");
fflush(gnu_pipe);
}
#if defined (_WIN32) && (_MSC_VER)
HWND hFrame = FindWindowEx(NULL,NULL,NULL,"XYPlot");
SwitchToThisWindow(hFrame, TRUE);
#endif
}
#endif//_MY_SIMPLE_PLOT_H_
- 来一个简单的hanning窗低通滤波器测试一下。
#include <fstream>
#include <math.h>
#include <vector>
#include "MySimplePlot.hpp"
typedef std::vector<double> DoubleVector;
const double PI = 3.141592653589793238462643383;
DoubleVector hanning_window(int N)
{
DoubleVector window(N);
int M = N - 1;
for(int i = 0; i < N; i++)
{
window[i] = 0.5 - 0.5*cos(2.0*PI*i/M);
}
return window;
}
DoubleVector lowpass_sinc(int N, double fc)
{
DoubleVector data(N);
int M = N - 1;
for(int i = 0; i < N; i++)
{
int n = i-M/2;
if(n==0) data[i] = 2.0*fc;
else data[i] = sin(n*2.0*PI*fc)/(n*PI);
}
return data;
}
DoubleVector convolution(const DoubleVector& x, const DoubleVector& h)
{
int x_length = x.size();
int h_length = h.size();
int length = x_length + h_length - 1;
DoubleVector result(length, 0.0);
for(int i = 0; i < length; i++)
{
for(int j = 0; j < h_length; j++)
{
if((i-j >= 0) && (i-j<x_length))
{
result[i] += h[j]*x[i-j];
}
}
}
return result;
}
DoubleVector inner_product(const DoubleVector& x, const DoubleVector& y)
{
int x_length = x.size();
int y_length = y.size();
if(x_length != y_length)
throw std::invalid_argument("inner_product() x.size() != y.size() ");
DoubleVector result(x_length);
for(int i = 0; i < x_length; i++)
{
result[i] = x[i]*y[i];
}
return result;
}
int main()
{
int data_length = 1000;
int filter_length = 155;
DoubleVector sine_before(data_length);
for(int i = 0; i <data_length; i++)
{
// 1V/1kHz + 0.1V/5kHz + 0.2V/10Hz
sine_before[i]=sin(2*PI*i*0.01) + 0.1*sin(2*PI*i*0.05) + 0.2*sin(2.0*PI*i*0.1);
}
DoubleVector win = hanning_window(filter_length);
DoubleVector sinc = lowpass_sinc(filter_length, 0.02);
DoubleVector h = inner_product(win, sinc);
DoubleVector sine_after = convolution(sine_before, h);
XYChart chart1;
chart1.title("FIR Filter");
chart1.add("before", data_length-filter_length+1, sine_before.data());
chart1.add("after", data_length-filter_length+1, sine_after.data());
chart1.plot();
XYChart chart2;
chart2.title("Filter Data Plot");
chart2.add("sinc", filter_length, sinc.data());
chart2.add("h", filter_length, h.data());
chart2.plot();
getchar();
return 0;
}
来两张生成的效果图: