3274|0

3836

帖子

19

TA的资源

纯净的硅(中级)

楼主
 

DSP28335通过FFT变换实现高频滤波 [复制链接]

本帖最后由 fish001 于 2019-11-5 21:17 编辑

由于AD采集到的高频开关噪声大,对数据处理存在极大干扰,经过示波器测量,发现噪声频率主要集中在80kHz,寻求老师帮助,老师推荐使用FFT滤波,有同学使用的是FIR滤波,通过Matlab中自带的FDA工具获得参数,效果在噪声幅值大于0.6V存在阶段性下滑,我使用的FFT滤波,由于网上资料有限,网上都是推荐的是直接使用FIR滤波这类的滤波器,FFT滤波可能会存在一些问题,但是目前对于本次项目,我还是觉得尝试一下FFT滤波,查看一下效果如何。
首先使用的Matlab进行仿真,代码如下

SampleFreq =    800000;            %为滤去80k高频 需大于2*80k的采样频率
S1_Freq    =    20000;            %模拟的初始信号频率
S2_Freq        =80000;        %高频开关噪声的频率
SAMPLE_NODES = 32;  
i = 1:SAMPLE_NODES;  

wt1=2*pi*i*S1_Freq;
wt1=wt1/SampleFreq;
wt2=2*pi*i*S2_Freq;
wt2=wt2/SampleFreq;

x = sin(wt1) +0.2*cos(wt2)+2; %获得20KHz基波+80KHz谐波信号,
                              %其中+2幅值是为了保证(实际采集到的)输入信号始终为正

subplot(2,2,1); plot(x);title('Inputs');  
axis([0 SAMPLE_NODES 0 5]);  


y = fft(x, SAMPLE_NODES);       %对x进行fft变换
subplot(2,2,2); plot(abs(y));title('FFT');  


axis([0 SAMPLE_NODES 0 80]);  

for m=3:31                                  %对fft变换后的x进行滤波
    y(m) = 0;
end

z = ifft(y, SAMPLE_NODES);  %对滤波后的y进行IFFT
subplot(2,2,3); plot(abs(z));title('IFFT');  
axis([0 SAMPLE_NODES 0 5]);  

subplot(2,2,4);; plot(abs(y));title('IFFT-FFT');  
axis([0 SAMPLE_NODES 0 80]);  

从而得到了以下FFT实现代码
头文件.h

/*
 * DSP_2833x_FFT.h
 *
 *  Created on: 2019年4月1日
 *      Author: ChenQingye
 */
/*     日期          实现功能
 *     0327    实现FFT变换,能成功滤掉高频部分,但还需完善:IFFT变换,2k频率以内进行滤波
 *             Graph坐标轴值:x轴:X(n) = (fs/N)*n        Y(n) = (N/2)*A    A为初始幅值默认1024
 *     0328    完善FFT功能,满足200kHZ采样频率采到80kHZ高频噪声。
 *     0401    自己编写ifft程序,需优化代码
 *     0402    使用256点的FFT时钟周期需310,048个时钟周期,总814,224时钟周期
 *             使用128点的FFT需要103,892个时钟周期,总386,111时钟周期
 *             使用64点的FFT需要69,980个时钟周期,总183,774时钟周期
 *             使用32点的FFT需要83,462个时钟周期,总83,462时钟周期
 *
 *
 *
 * */
#ifndef DSP_2833X_FFT_H_
#define DSP_2833X_FFT_H_


#include "DSP28x_Project.h"     // Device Headerfile and Examples Include File
#include"math.h"
#define SampleFreq     1000000            // 为滤去80k高频 需大于2*80k的采样频率
#define S1_Freq        20000            // 模拟的初始信号频率
#define S2_Freq        80000            // 高频开关噪声的频率
#define PI 3.1415926
#define SAMPLENUMBER 32            // 0-256个FFT样本点涵盖了一个低频信号0-FS的所有频率信息,FS是采样频率。
#define M log(SAMPLENUMBER)/log(2)     // 变更FFT点数修改fft子程序方便的需要    2^8 = 256->M = 8;log(SAMPLENUMBER)/log(2)

void InitForFFT();
void FFT256(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER]);//新加x7
void FFT128(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER]);
void FFT64(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER]);//去掉x6
void FFT32(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER]);//去掉x5
void IFFT();
void FilterFFT();//实现FFT滤波器
void MakeFFTOut();
void MakeIFFTOut();
void CLearFFtData();
//void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]);

extern int16 INPUT[SAMPLENUMBER],FFTDATA[SAMPLENUMBER],IFFTDATA[SAMPLENUMBER];
extern float32 DATAR[SAMPLENUMBER],DATAI[SAMPLENUMBER];
extern float32 sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];

#endif /* DSP_2833X_FFT_H_ */


实现代码

/*
 * DSP_2833x_FFT.c
 *
 *  Created on: 2019年4月3日
 *      Author: ChenQingye
 */
#include "DSP_2833x_FFT.h"


int16 INPUT[SAMPLENUMBER],FFTDATA[SAMPLENUMBER],IFFTDATA[SAMPLENUMBER];
float32 DATAR[SAMPLENUMBER],DATAI[SAMPLENUMBER];
float32 sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];


void InitForFFT()//蝶形运算系数表计算
{
    int16 i;

    float32 wt1,wt2;


    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        // 生成蝶形运算需要的三角表
        sin_tab=sin(PI*2*i/SAMPLENUMBER);
        cos_tab=cos(PI*2*i/SAMPLENUMBER);        // f = w/2pi =1hz

        // 生成模拟信号
        wt1=2*PI*i*S1_Freq;
        wt1=wt1/SampleFreq;

        wt2=2*PI*i*S2_Freq;
        wt2=wt2/SampleFreq;

        INPUT=sin(wt1)*1024+0.2*cos(wt2)*1024+1024;        //f = w/2pi = 3hz
    }
    //初始化输入数据以及清空输出数组
    for ( i=0;i<SAMPLENUMBER;i++ )
        {
            DATAR = INPUT;
            DATAI = 0.0f;
            FFTDATA = 0.0f;
            IFFTDATA = 0.0f;
        }
}

void FilterFFT()//实现FFT滤波器
{
    int16 i = 0;
    int16 upLimit,downLimit;
    switch(SAMPLENUMBER)
        {
        case 32:
            upLimit = 31;
            downLimit = 3;
            break;
        case 64:
            upLimit = 62;
            downLimit = 4;
            break;
        case 128:
            upLimit = 120;
            downLimit = 9;
            break;
        case 256:
            upLimit = 240;
            downLimit = 12;
            break;
        default:
            break;
        }

    for(i = downLimit;i < upLimit;i++)//100kHZ采样 256点
    {
        DATAR = 0;
        DATAI = 0;
    }
}
void IFFT()
{
    int16 k = 0;

    for (k=0; k<=SAMPLENUMBER-1; k++) {
        DATAI[k] = -DATAI[k];
    }

    switch(SAMPLENUMBER)
    {
    case 32:
        FFT32(DATAR,DATAI);
        break;
    case 64:
        FFT64(DATAR,DATAI);
        break;
    case 128:
        FFT128(DATAR,DATAI);
        break;
    case 256:
        FFT256(DATAR,DATAI);
        break;
    default:
        break;
    }
//
//    if(SAMPLENUMBER == 256)/* using FFT */
//        FFT256(DATAR,DATAI);
//    else if(SAMPLENUMBER == 128)
//        FFT128(DATAR,DATAI);
//    else if(SAMPLENUMBER == 64)
//        FFT64(DATAR,DATAI);
//    else if(SAMPLENUMBER == 32)
//            FFT32(DATAR,DATAI);

    for (k=0; k<=SAMPLENUMBER-1; k++) {
        DATAR[k] = DATAR[k] / SAMPLENUMBER;
        DATAI[k] = -DATAI[k] / SAMPLENUMBER;
    }

}
void MakeFFTOut()
{
    Uint16 i;
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        FFTDATA=sqrt(DATAR*DATAR+DATAI*DATAI);
    }
}

void MakeIFFTOut()
{
    Uint16 i;
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        IFFTDATA=sqrt(DATAR*DATAR+DATAI*DATAI);
    }
}


void FFT256(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER])//新加x7    310,048个时钟周期
{
    int16 x0,x1,x2,x3,x4,x5,x6,xx,x7;//x7
    int16 i,j,k,b,p,L;
    float32 TR,TI,temp;

    /********** 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;x7=(i/128)&0x01;//x7
        xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7;
        dataI[xx]=dataR;
    }
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        dataR=dataI; dataI=0;
    }

    /************** following code FFT *******************/
    for ( L=1;L<=M;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=M-L;
            while ( i>0 ) /* p=pow(2,M-L)*j; */
            {
                p=p*2; i--;
            }
            p=p*j;
            for ( k=j;k<SAMPLENUMBER;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;i++ )
    {
        DATAR = dataR;
        DATAI = dataI;
    }
} /* END FFT */

void FFT128(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER])    //103,892个时钟周期
{
    int16 x0,x1,x2,x3,x4,x5,x6,xx;
    int16 i,j,k,b,p,L;
    float32 TR,TI,temp;

    /********** 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;
    }
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        dataR=dataI; dataI=0;
    }

    /************** following code FFT *******************/
    for ( L=1;L<=M;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=M-L;
            while ( i>0 ) /* p=pow(2,7-L)*j; */
            {
                p=p*2; i--;
            }
            p=p*j;
            for ( k=j;k<SAMPLENUMBER;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;i++ )
    {
        DATAR = dataR;
        DATAI = dataI;
    }
} /* END FFT */

void FFT64(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER])    //103,892个时钟周期
{
    int16 x0,x1,x2,x3,x4,x5,xx;
    int16 i,j,k,b,p,L;
    float32 TR,TI,temp;

    /********** following code invert sequence ************/
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        x0=x1=x2=x3=x4=x5=0;
        x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01;
        xx=x0*32+x1*16+x2*8+x3*4+x4*2+x5;
        dataI[xx]=dataR;
    }
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        dataR=dataI; dataI=0;
    }

    /************** following code FFT *******************/
    for ( L=1;L<=M;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=M-L;
            while ( i>0 ) /* p=pow(2,7-L)*j; */
            {
                p=p*2; i--;
            }
            p=p*j;
            for ( k=j;k<SAMPLENUMBER;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;i++ )
    {
        DATAR = dataR;
        DATAI = dataI;
    }
} /* END FFT */

void FFT32(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER])    //103,892个时钟周期
{
    int16 x0,x1,x2,x3,x4,xx;
    int16 i,j,k,b,p,L;
    float32 TR,TI,temp;

    /********** following code invert sequence ************/
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        x0=x1=x2=x3=x4=0;
        x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01;
        xx=x0*16+x1*8+x2*4+x3*2+x4;
        dataI[xx]=dataR;
    }
    for ( i=0;i<SAMPLENUMBER;i++ )
    {
        dataR=dataI; dataI=0;
    }

    /************** following code FFT *******************/
    for ( L=1;L<=M;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=M-L;
            while ( i>0 ) /* p=pow(2,7-L)*j; */
            {
                p=p*2; i--;
            }
            p=p*j;
            for ( k=j;k<SAMPLENUMBER;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;i++ )
    {
        DATAR = dataR;
        DATAI = dataI;
    }
} /* END FFT */


最终采用的是800KHz采样频率,32个点。由于工程中PID调节的速度有限,32个点能满足滤波要求。理论中点数越多,滤波效果越好,800KHz/32点=25kHz,表示25KHz以下的频率无法滤除,而1.6MHz中50KHz以下都无法滤除,才导致滤波效果不好
 

点赞 关注
 

回复
举报
您需要登录后才可以回帖 登录 | 注册

随便看看
查找数据手册?

EEWorld Datasheet 技术支持

相关文章 更多>>
关闭
站长推荐上一条 1/9 下一条

 
EEWorld订阅号

 
EEWorld服务号

 
汽车开发圈

About Us 关于我们 客户服务 联系方式 器件索引 网站地图 最新更新 手机版

站点相关: 国产芯 安防电子 汽车电子 手机便携 工业控制 家用电子 医疗电子 测试测量 网络通信 物联网

北京市海淀区中关村大街18号B座15层1530室 电话:(010)82350740 邮编:100190

电子工程世界版权所有 京B2-20211791 京ICP备10001474号-1 电信业务审批[2006]字第258号函 京公网安备 11010802033920号 Copyright © 2005-2025 EEWORLD.com.cn, Inc. All rights reserved
快速回复 返回顶部 返回列表