/****************************************************************************
*
*文件名：AD02.c
*
*功 能：实现采样数据的快速傅立叶变换，并进行功率谱计算，并对频谱进行修正
*
*作 者:
*
****************************************************************************/
/*
说 明：
输 入：
Sampling:采集点数
fs:采样频率
fWaveR:采集数据实部
fWaveI:采集数据虚部
//W:修正后的幅值
//power:功率谱
//num:谱线数
输 出:
PinPu:频谱计算后幅值
power:功率谱
NengLiang0:0~300Hz能量
NengLiang1:300~500Hz能量
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926
#define SAMPLENUMBER 1024//采样点数
#define FIRNUMBER 16//FIR窗体长度
//FFT的参数变量
float fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER];
float w[SAMPLENUMBER/2];//频谱幅值
//float W[SAMPLENUMBER/2];//修正后频谱幅值
float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];//
/************功率谱计算************************/
//float max; //功率谱的最大值
float power[SAMPLENUMBER/2];//功率谱
//float powerf[SAMPLENUMBER/2];//功率谱峰值谱线号
//int num=0;//峰值个数
//文本数据显示内容
float BoXing[SAMPLENUMBER];
float PinPu[SAMPLENUMBER/2];
float PinPu50,PinPu100,PinPu150,PinPu300;
float NengLiang0,NengLiang1;
//FIR的参数
float fHn[FIRNUMBER]=
{
5.882744631e-019f,-0.0001204110522f, 0.004326796625f, 0.004784306977f, -0.03239288554f,
-0.04123774916f, 0.1447458863f, 0.4198940396f, 0.4198940396f, 0.1447458863f,
-0.04123774916f, -0.03239288554f, 0.004784306977f, 0.004326796625f,-0.0001204110522f,
5.882744631e-019f
};
float fXn[FIRNUMBER]={ 0.0 };
float fIn[SAMPLENUMBER],fOut[SAMPLENUMBER];
int nIn,nOut;
void MakeWave()
{
int i;
for ( i=0;i<SAMPLENUMBER;i++ )
{
fWaveR[i]=float(1+sin(PI*2*i*0.002/5.12)+sin(PI*2*i*0.1/5.12)+2*sin(PI*2*i*1.5/5.12));//sin函数为弧度；
fIn[i]=BoXing[i]=fWaveR[i];
fWaveI[i]=0;
//ShiJian[i]=(float)i*8000/SAMPLENUMBER;
//shuju.Wave[i]=float(1+sin(PI*4*i/SAMPLENUMBER)+sin(PI*8*i/SAMPLENUMBER));//sin函数为弧度
}
}
/*****快速傅里叶变换初始化*************/
void InitForFFT()
{
int i;
for ( i=0;i<SAMPLENUMBER;i++ )
{
sin_tab[i]=(float)sin(PI*2*i/SAMPLENUMBER);
cos_tab[i]=(float)cos(PI*2*i/SAMPLENUMBER);
}
}
/*****快速傅里叶变换函数************/
void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER])
{
int x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,xx;
int i,j,k,b,p,L;
float TR,TI,temp;
/********** following code invert sequence ************/
for ( i=0;i<SAMPLENUMBER;i++ )
{
x0=x1=x2=x3=x4=x5=x6=x7=x8=x9=0;
x0=i&0x001; x1=(i/2)&0x0001; x2=(i/4)&0x001; x3=(i/8)&0x001;x4=(i/16)&0x001; x5=(i/32)&0x001; x6=(i/64)&0x001;x7=(i/128)&0x001;x8=(i/256)&0x001;x9=(i/512)&0x001;
xx=x0*512+x1*256+x2*128+x3*64+x4*32+x5*16+x6*8+x7*4+x8*2+x9;
dataI[xx]=dataR[i];
}
for ( i=0;i<SAMPLENUMBER;i++ )
{
dataR[i]=dataI[i]; dataI[i]=0;
}
/************** following code FFT *******************/
for ( L=1;L<=10;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=10-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/2;i++ )
{
if(i==0)
w[i]=(float)sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i])/SAMPLENUMBER;//直流频谱幅值
else
w[i]=(float)sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i])*2/SAMPLENUMBER;//频谱幅值
PinPu[i]=w[i];
}
/************功率谱计算************************/
for(i=0;i<SAMPLENUMBER/2;i++)
{
power[i]=w[i]*w[i];
}
/************能量计算************************/
NengLiang0=0;
NengLiang1=0;
//0~300Hz能量
for(i=0;i<60;i++)
{
NengLiang0=NengLiang0+w[i];
}
//300~500Hz能量
for(i=60;i<100;i++)
{
NengLiang1=NengLiang1+w[i];
}
} /* END FFT */
//对故障进行分析
void GuZhangFenXi()
{
if((PinPu[20]/PinPu100>1.4)||(4.07>NengLiang0/NengLiang1>3.67))
{
printf("匝间短路故障\n");
}
if(((PinPu[10]/PinPu50>8)&&(PinPu[30]/PinPu150>18))||(19.69>NengLiang0/NengLiang1>19.29))
{
printf("绕组变形故障\n");
}
if((PinPu[60]/PinPu300>1.4)||(2.28>NengLiang0/NengLiang1>1.88))
{
printf("铁心松动故障\n");
}
}
//FIR输入函数
float InputWave()
{
int i=0;
for ( i=FIRNUMBER-1;i>0;i-- )
{
fXn[i]=fXn[i-1];
}
fXn[0]=fIn[nIn];
nIn++;
return(fXn[0]);
}
//FIR函数
float FIR()
{
int i=0;
float fSum;
fSum=0;
for ( i=0;i<FIRNUMBER;i++ )
{
fSum+=(fXn[i]*fHn[i]);
}
return(fSum);
}
void main()
{
int i;
FILE * fp;
nIn=0;
nOut=0;
//生成波形
MakeWave();
//滤波前的FFT分析
InitForFFT();
//FFT计算
FFT(fWaveR,fWaveI); //@@@@@@@@@@@@@@@@@@@2
GuZhangFenXi();
//采样数据的存储
if((fp=fopen("F:\\波形.txt","w"))==NULL)
{
printf("文件不能打开!\n");
exit(0);
} //为输出打开一个文本文件
for( i = 0; i <SAMPLENUMBER; ++i )
{
fprintf(fp,"%f\n",BoXing[i]);
}
fclose(fp);
//FFT分析数据的存储
if((fp=fopen("F:\\频谱1.txt","w"))==NULL)
{
printf("文件不能打开!\n");
exit(0);
} //为输出打开一个文本文件
for( i = 0; i <SAMPLENUMBER/2; ++i )
{
fprintf(fp,"%f\n",PinPu[i]);
}
fprintf(fp,"%f\n",PinPu[300]);
fclose(fp);
//FIR滤波
for(i=0;i<SAMPLENUMBER;i++)
{
InputWave();
fOut[nOut]=FIR();
nOut++;
}
//滤波后的FFT分析
FFT(fOut,fWaveI); //@@@@@@@@@@@@@@@@@@@2
//FFT分析数据的存储
GuZhangFenXi();
if((fp=fopen("F:\\频谱2.txt","w"))==NULL)
{
printf("文件不能打开!\n");
exit(0);
} //为输出打开一个文本文件
for( i = 0; i <SAMPLENUMBER/2; ++i )
{
fprintf(fp,"%f\n",PinPu[i]);
}
fprintf(fp,"%f\n",PinPu[300]);
fclose(fp);
}
