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;
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]);  

32点800Khz采样Matlab仿真图
在经过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;//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[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库进行运算,会增添许多不必要的运算时间,这点折腾了我好一会儿。

  • 2
    点赞
  • 39
    收藏
    觉得还不错? 一键收藏
  • 7
    评论
DSP28335是一款数字信号处理器,可以进行FFT(快速傅里叶变换)操作。引用\[1\]中提到,虽然例程中没有提供FFT的例程,但是可以使用C2000 ware库中的C2000系列库来进行FFT操作。同时,引用\[2\]中介绍了DSP28335的C28X_FPU_LIB.lib库中的CFFT_f32函数,该函数用于计算复数域的FFT。该函数需要传入一个CFFT_F32_STRUCT结构体作为参数,该结构体中包含了傅里叶变换的阶数和转化因子的计算。引用\[3\]中提到,可以使用CFFT_f32_sincostable函数来计算傅里叶变换的转化因子,只需要在结构体中对stage和FFTSize进行赋值即可。因此,通过使用C2000 ware库和C28X_FPU_LIB.lib库中的函数,可以在DSP28335上进行FFT变换。 #### 引用[.reference_title] - *1* [DSP28335FFT傅里叶变换](https://blog.csdn.net/weixin_30390075/article/details/95330151)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^koosearch_v1,239^v3^insert_chatgpt"}} ] [.reference_item] - *2* *3* [DSP:关于28335FFT](https://blog.csdn.net/Jason_nuc/article/details/81975630)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^koosearch_v1,239^v3^insert_chatgpt"}} ] [.reference_item] [ .reference_list ]

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值