DSP28335通过FFT变换实现高频滤波

这次项目中，由于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;
endz = 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]);

在经过matlab仿真后,发现能达到预期滤波效果，动手！
在网上查阅各种资料，发现本科期间学过fft但是没怎么使用过，对FFT原理不太有印象了，恰好发现一篇神文：
https://zhuanlan.zhihu.com/p/19763358
对知识点进行了巩固哈哈哈哈

从而得到了以下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[i]=sin(PI*2*i/SAMPLENUMBER);cos_tab[i]=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[i]=sin(wt1)*1024+0.2*cos(wt2)*1024+1024;       //f = w/2pi = 3hz}//初始化输入数据以及清空输出数组for ( i=0;i<SAMPLENUMBER;i++ ){DATAR[i] = INPUT[i];DATAI[i] = 0.0f;FFTDATA[i] = 0.0f;IFFTDATA[i] = 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[i] = 0;DATAI[i] = 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[i]=sqrt(DATAR[i]*DATAR[i]+DATAI[i]*DATAI[i]);}
}void MakeIFFTOut()
{Uint16 i;for ( i=0;i<SAMPLENUMBER;i++ ){IFFTDATA[i]=sqrt(DATAR[i]*DATAR[i]+DATAI[i]*DATAI[i]);}
}void FFT256(float32 dataR[SAMPLENUMBER],float32 dataI[SAMPLENUMBER])//新加x7 310,048个时钟周期
{int16 x0,x1,x2,x3,x4,x5,x6,xx,x7;//x7int16 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;//x7xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7;dataI[xx]=dataR[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=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[i] = dataR[i];DATAI[i] = dataI[i];}
} /* 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[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=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[i] = dataR[i];DATAI[i] = dataI[i];}
} /* 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[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=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[i] = dataR[i];DATAI[i] = dataI[i];}
} /* 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[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=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[i] = dataR[i];DATAI[i] = dataI[i];}
} /* END FFT */

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

参考文章：
https://blog.csdn.net/wang_0712/article/details/80543801
http://bbs.21ic.com/icview-818068-1-1.html
https://blog.csdn.net/dreamandxiaochouyu/article/details/45392723
以上。
另：在优化代码时，确定好点数后，将M的值log(SAMPLENUMBER)/log(2)修改为确定的一个数而不是调用math库进行运算，会增添许多不必要的运算时间，这点折腾了我好一会儿。

DSP28335通过FFT变换实现高频滤波相关推荐

干货 | 使用FFT变换自动去除图像中严重的网纹
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达本文转自:opencv学堂这个课题在很久以前就已经有所接触,不过 ...
【算法随记五】使用FFT变换自动去除图像中严重的网纹。
这个课题在很久以前就已经有所接触,不过一直没有用代码去实现过.最近买了一本<机器视觉算法与应用第二版>书,书中再次提到该方法:使用傅里叶变换进行滤波处理的真正好处是可以通过使用定制的滤波器 ...
数字图像处理100问—33 傅立叶变换——低通滤波
提示:内容整理自:https://github.com/gzr2017/ImageProcessing100Wen CV小白从0开始学数字图像处理 33 傅立叶变换--低通滤波将imori.jpg灰 ...
功率谱和FFT变换频谱
功率谱: nfft=length(total_wave); window1=hamming(100); %海明窗 noverlap=20; %数据无重叠 range='onesided'; %频率间隔 ...
使用FFT变换自动去除图像中严重的网纹
本篇博文来自博主Imageshop,打赏或想要查阅更多内容可以移步至Imageshop. 转载自:https://www.cnblogs.com/Imageshop/p/11625895.html侵删 ...
语音信号分析之FFT变换过程
语音信号分析从FFT变换开始概述语音信号是人为经过声学设备采集转换而来的便于编码.储存以及做声学处理用途抽象而成的信号.首先一段WAV音频拥有固定的几个参数用来规定语音信号的属性.正确的区分通道数 ...
深入浅出解释FFT（六）——深入理解fft变换
(如需交流,请关注公众号:神马观止) FFT(FastFourier Transform,快速傅立叶变换)是离散傅立叶变换的快速算法,也是我们在数字信号处理技术中经常会提到的一个概念.在大学的理工科课 ...
OpenCV使用Laplacian filtering和距离变换以及Laplacian滤波对重叠对象进行分段的实例(附完整代码)
OpenCV使用Laplacian filtering和距离变换以及Laplacian滤波对重叠对象进行分段的实例 OpenCV使用Laplacian filtering和距离变换以及Laplacia ...
电路中滤波电容和退耦电容_详解电源滤波电路中的高频滤波电容电路
图2-12所示是电源滤波电路中的高频滤波电路.电路中,一个容量很大的电解电容C1(2200µF)与一个容量很小的电容C2(0.01µF)并联,C2是高频滤波电容,用来进行高频成分的滤波,这种一大一小两 ...

DSP28335通过FFT变换实现高频滤波

DSP28335通过FFT变换实现高频滤波相关推荐

最新文章

热门文章