AD采样实现AC计量之Matlab 绘制C函数图像篇（一）

笔者按：

        十一长假来了，感觉无聊呀，就看一点FFT的资料，研究一下，其实有时候在学习的时候你也能得到一种乐趣，特别是结果出来的时候，那时候会有一种成就感，以前在做智能小车的时候要调试PID当时想过用matlab仿真来计算pid的各种参数，但限于当时时间和技术水平有限最后没有能搭建起这个平台，但现在公司有一个项目的一部分是dsp采集交流电计算出频率、幅值、有效值、相角、谐波等，作者一直不情愿搞dsp这一块，arm才是笔者的钟爱，不过这个项目牵扯到dsp的一些算法，复杂的算法也是dsp的有优势和玩头的地方，笔者的前一份工作是做电表的，对这些东西虽然不陌生但是以前的计量参数，全是由计量芯片实现的，计量芯片内部也是dsp和高精度ad实现的，这个东西做好了对笔者也是极大的提高，于是笔者决定利用假期的时间，一点点把这个硬柿子拿下，年轻的时候不要怕，有什么新的事物都可以尝试，碰到自己没有头绪的东西，要有坚定的信念，很多东西就是赢在坚定的信念上的，这些技术上的问题都不是核心问题，不懂可以学，但一个整体的格局要上去，一个人的格局很多时候决定了你的最高点，笔者准备把这些东西记录下来，一方面使后来人少走弯路，一方面留下一点青春的印迹。

我们来看一下如何利用matlab绘制出c函数的图像便于一些算法的调试：

在matlab中MyFFT.c内容如下：

#include "mex.h"    //头文件必须包含mex.h#include "math.h"#define PI 3.1415926#define SAMPLENUMBER 128double FFT(double w[],double dataR[],double n);//c语言到matlab变换，以mexFunction命名void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[]){    double *pw;    double *pda;        plhs[0]=mxCreateDoubleMatrix(1,SAMPLENUMBER,mxREAL);    pw=mxGetPr(plhs[0]);    pda=mxGetPr(prhs[0]);    FFT(pw,pda,0);}//C语言函数  double FFT(double w[],double dataR[],double n){    int x0, x1, x2, x3, x4, x5, x6, xx;    int i, j, k, b, p, L;    double TR, TI, temp;    double sin_tab[SAMPLENUMBER], cos_tab[SAMPLENUMBER];    double dataI[SAMPLENUMBER];    for (i = 0; i < SAMPLENUMBER; i++)    {        sin_tab[i] = sin(PI *2 * i / SAMPLENUMBER);        cos_tab[i] = cos(PI *2 * i / SAMPLENUMBER);    }    /********** following code invert sequence ************/    for (i = 0; i < SAMPLENUMBER; i++)    {        x0 = x1 = x2 = x3 = x4 = x5 = x6 = 0;        x0 = i &0x01;        x1 = (i / 2) &0x01;        x2 = (i / 4) &0x01;        x3 = (i / 8) &0x01;        x4 = (i / 16) &0x01;        x5 = (i / 32) &0x01;        x6 = (i / 64) &0x01;        xx = x0 * 64+x1 * 32+x2 * 16+x3 * 8+x4 * 4+x5 * 2+x6;        dataI[xx] = dataR[i];    }    for (i = 0; i < SAMPLENUMBER; i++)    {        dataR[i] = dataI[i];        dataI[i] = 0;    }    /************** following code FFT *******************/    for (L = 1; L <= 7; L++)    {         /* for(1) */        b = 1;        i = L - 1;        while (i > 0)        {            b = b * 2;            i--;        } /* b= 2^(L-1) */        for (j = 0; j <= b - 1; j++)         /* for (2) */        {            p = 1;            i = 7-L;            while (i > 0)             /* p=pow(2,7-L)*j; */            {                p = p * 2;                i--;            }            p = p * j;            for (k = j; k < 128; k = k + 2 * b)             /* for (3) */            {                TR = dataR[k];                TI = dataI[k];                temp = dataR[k + b];                dataR[k] = dataR[k] + dataR[k + b] *cos_tab[p] + dataI[k + b] *sin_tab[p];                dataI[k] = dataI[k] - dataR[k + b] *sin_tab[p] + dataI[k + b] *cos_tab[p];                dataR[k + b] = TR - dataR[k + b] *cos_tab[p] - dataI[k + b] *sin_tab[p];                dataI[k + b] = TI + temp * sin_tab[p] - dataI[k + b] *cos_tab[p];            } /* END for (3) */        } /* END for (2) */    } /* END for (1) */    for (i = 0; i < SAMPLENUMBER / 2; i++)    {        w[i] = sqrt(dataR[i] *dataR[i] + dataI[i] *dataI[i]);    }} /* END FFT */

上面的FFT函数是Ti的基2离散FFT算法，memFunction是matlab提供的转化接口函数，关于用法笔者之前的博客讲到过，不清楚的请补课。

在命令窗运行 

>>mex MyFFT.c

编写lastT.m文件，内容如下：

fs=200; %设定采样频率N=128;n=0:N-1;t=n/fs;f0=50; %正弦频率为10HZ% 正弦信号发生x=(sin(2*pi*f0*t));figure(1);subplot(231); %第一个子图plot(t,x); %作正弦信号的时域波形xlabel('t');ylabel('y');title('y=2*pi*10t时域');grid;% 进行FFT变换并做频谱图x1=(sin(2*pi*f0*t));%x1=square(2*pi*f0*t,50);y1=MyFFT(x1,N); %进行自己的fft转换%y1=fft(x1,N);mag1=abs(y1)*2/N; %幅值m1=length(y1); f1=(0:m1/2-1)'*fs/m1; %进行对应的频率转换%x2=(sin(2*pi*f0*t));% x2=square(2*pi*f0*t,50);% y2=fft(x2,N); %标准fft% mag2=abs(y2); % m2=length(y); % f2=(0:m2/2-1)'*fs/m2; %进行对应的频率转换figure(1);subplot(232);%plot(f1,mag1(1:m1/2),'g',f2,mag2(1:m2/2),'-.r');%做my频谱图plot(f1,mag1(1:m1/2),'g');%做频谱图%plot(f2,mag2(1:m/2),'r');%做频谱图axis([0,100,0,1.1]);xlabel('频率(Hz)');ylabel('幅值');title('y=2*pi*10t频域');grid;

运行结果图

经过fft转化的会有个幅度的数组，利用这个数组可以换算出真实的幅度值和频率值

假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。

其实这样出来的精度偏差还是很大的，第一是受采样点数的限制，第二采样的频率是真实频率2的整数倍的时候幅值的偏差相对较小。