本帖最后由 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以下都无法滤除,才导致滤波效果不好
|