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 128


double 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的整数倍的时候幅值的偏差相对较小。