📄 fft.cpp
字号:
#include "stdafx.h"
#include "Fft.h"
#include "math.h"
const double MOD_MAX = 65535.0;
const double DISP_MAX = 1.0/255.0;
extern FILE *fp;
//FFT运算必须参数
int fft_point,fft_order,fft_divide,fft_window,fft_scale;
bool fft_cover;
float filter[7];//FIR滤波参数
extern FILE *fpIandQ;
extern bool m_bIqWrite;
double prFilter[256],piFilter[256];
/**************************************************************************************
0 RectangleWindow矩形窗
FFT变换结果为对称型，矩形窗是使信号突然截断，旁瓣会很大，且衰减较慢，旁瓣的第一个负
峰值为主瓣的21％，第一个正峰值为主瓣的12.6％，第二个负峰值为主瓣的9％，效果一般，泄
漏较大。
**************************************************************************************/
double WINAPI RectangleWindow(int t)
{
double wt;
wt=1.0;
return wt;
}
/**************************************************************************************
1 TriangleWindow三角窗,也称费杰（Fejer）窗，Bartlett
**************************************************************************************/
double WINAPI TriangleWindow(int t)
{
double wt;
wt=1-t/fft_point;
return wt;
}
/**************************************************************************************
2 HanningWindwo汉宁窗，即升余弦窗
**************************************************************************************/
double WINAPI HanningWindow(int t)
{
double wt;
wt=(1-cos(2*PI*t/fft_point))/2;
return wt;
}
/**************************************************************************************
3 HammingWindow海明窗，即改进的升余弦窗
**************************************************************************************/
double WINAPI HammingWindow(int t)
{
double wt;
wt=0.54-0.46*cos(2*PI*t/fft_point);
return wt;
}
/**************************************************************************************
4 BlackmanWindow布来克曼窗，即三阶升余弦窗
**************************************************************************************/
double WINAPI BlackmanWindow(int t)
{
double wt;
wt=0.42-0.5*cos(2*PI*t/fft_point)+0.08*cos(4*PI*t/fft_point);
return wt;
}
/**************************************************************************************
5 CosgradeWindow余弦坡度窗
**************************************************************************************/
double WINAPI CosgradeWindow(int t)
{
double wt;
if(t<= int(4*fft_point/5) )
wt=1.0;
else
wt=(1+cos(5*PI*t/fft_point))/2;
return wt;
}
/**************************************************************************************
6 ParzenWindow帕曾窗
**************************************************************************************/
double WINAPI ParzenWindow(int t)
{
double wt;
if(t<= int(fft_point/2) )
wt=1-6*pow(t/fft_point,2)+6*pow(t/fft_point,3);
else
wt=2*pow(1-t/fft_point,3);
return wt;
}
/**************************************************************************************
7 ExponentWindow指数窗
**************************************************************************************/
double WINAPI ExponentWindow(int t)
{
double wt,alf;
alf=0.0001;
wt=exp((-1.0)*alf*t);
return wt;
}
/**************************************************************************************
8 GaussWindow高斯窗
**************************************************************************************/
double WINAPI GaussWindow(int t)
{
double wt,alf;
alf=0.0001;
wt=exp(-1*alf*t*t);
return wt;
}
/**************************************************************************************
9 KaiserWindow凯泽窗,
凯泽在1966（1974）发现，利用第一类零阶修正（变形）贝赛尔函数可以构成一种近似最佳的窗
最主要参数：beat－可同时调整主瓣宽度和旁瓣，beat越大，窗越窄，频谱旁瓣越小，但主瓣相应增加
beta=0相当与矩形窗
**************************************************************************************/
double WINAPI KaiserWindow(int t)
{
double wt,alfa,beta;
alfa=(fft_point-1)/2;
beta=0;
//Bessel零阶第一类贝塞尔函数
wt=Jim_Bessel0_R(0,beta*sqrt(1-pow((t-alfa)/alfa,2)))/Jim_Bessel0_R(0,beta);
return wt;
}
/**************************************************************************************
10 ChebWindow契比雪夫窗
**************************************************************************************/
double WINAPI ChebWindow(int t)
{
double wt;
wt=1.0;
return wt;
}
/**************************************************************************************
11 BartlettWindow巴特利特窗
**************************************************************************************/
double WINAPI BartlettWindow(int t)
{
double wt;
wt=1.0;
return wt;
}
bool WINAPI fftInit(int point,int order,int divide,bool cover,int scale)
{
fft_point=point;
fft_order=order;
fft_divide=divide;
fft_cover=cover;
fft_scale=scale;
//采用8阶的FIR滤波
Jim_FirFilter(7,3,50,7800/2.0,7800.0,1,filter);
return TRUE;
}
BOOL WINAPI fftEnd()
{
return TRUE;
}
//傅立叶变换
BOOL WINAPI fftTransform(VOID *inBuf,VOID *outBuf,int windows,int scale,int nProbe,int nWork)
{
BYTE *srcPtr = (BYTE*)inBuf;//数据源1024字节
BYTE *dstPtr = (BYTE*)outBuf;//数据果256字节
//double *testPtr =(double *)testBuf;
WORD Ivalue,Qvalue;
//WORD IvalueRev,QvalueRev;
unsigned char flagQ,flagI;
//double alfa,e;
int i,j,Iorg=0,Qorg=0;
double mod = 0;
if(nProbe==LEFT_PROBE)
{
flagQ=0x30,flagI=0x20;
}
else
{
flagQ=0x10,flagI=0x0;
}
//每次取4个字节数据，分离I/Q分量，判断db(12)，0表示I分量，1表示Q分量
if (((srcPtr[1]&0xf0) ==flagQ) && ((srcPtr[3]&0xf0) ==flagI))
{
Iorg = 2;
Qorg = 0;
}
if (((srcPtr[1]&0xf0) ==flagI) && ((srcPtr[3]&0xf0) ==flagQ))
{
Iorg = 0;
Qorg = 2;
}
double pr[256],pi[256],w;
//计算自相关函数和互相关函数
//double Rii=0.0,Rqq=0.0,Riq=0.0;
for ( i=0;i<fft_point;i++ )
{
//更改存储的偏移地址分离I/Q，
Ivalue = *((short*)(srcPtr + i*4 + Iorg));
Qvalue = *((short*)(srcPtr + i*4 + Qorg));
//~位非运算，低12位位有效数据
pr[i]=double(Ivalue & 0xfff);
pi[i]=double(Qvalue & 0xfff);
/*Rii=Rii+double(pr[i]*pr[i]);
Rqq=Rqq+double(pi[i]*pi[i]);
Riq=Riq+double(pr[i]*pi[i]);*/
}
/*Rii=Rii/fft_point;
Rqq=Rqq/fft_point;
Riq=Riq/fft_point;
alfa=asin(Riq/sqrt(Rii*Rqq));
e=sqrt(Rqq/Rii)-1;*/
for ( i=0;i<fft_point;i++ )
{
//防止谱泄漏，进行加窗处理
switch(windows)
{
case WND_RECTANGLE: w=RectangleWindow(i); break;
case WND_TRIANGLE: w=TriangleWindow(i); break;
case WND_HANNING: w=HanningWindow(i); break;
case WND_HAMMING: w=HammingWindow(i); break;
case WND_BLACKMAN: w=BlackmanWindow(i); break;
case WND_COSGRADE: w=CosgradeWindow(i); break;
case WND_PARZEN: w=ParzenWindow(i); break;
case WND_EXPONENT: w=ExponentWindow(i); break;
case WND_GAUSS: w=GaussWindow(i); break;
case WND_KAISER: w=KaiserWindow(i); break;
case WND_CHEB: w=ChebWindow(i); break;
case WND_BARTLETT: w=BartlettWindow(i); break;
default: w=RectangleWindow(i); break;
}
/*pr[i]=((1+e)*cos(alfa)*pr[i])/((1+e)*cos(alfa))*w;
pi[i]=(pi[i]-(1+e)*sin(alfa)*pr[i])/((1+e)*cos(alfa))*w;*/
pr[i]=pr[i]*w;
pi[i]=pi[i]*w;
prFilter[i]=piFilter[i]=0.0;
for(j=0;j<7;j++)
{
if(i+j>255)
{
prFilter[i]=prFilter[i];
piFilter[i]=piFilter[i];
}
else
{
prFilter[i]=prFilter[i]+pr[i+j]*w*filter[j];
piFilter[i]=piFilter[i]+pi[i+j]*w*filter[j];
}
}
//fprintf(fp,"%d,%f\n",i,filter[i]);
//pr[i]=pr[i]*filter[i];
//pi[i]=pi[i]*filter[i];
}
//fclose(fpIandQ);
//fft处理
//Jim_FFT(pr,pi,fft_point,fft_order,0,1);
Jim_FFT(prFilter,piFilter,fft_point,fft_order,0,1);
/*for(int t=0;t<256;t++)
{
fprintf(fpIandQ,"%d\n",int(prFilter[t]));
}*/
//归一化
//Jim_Unitary(pr,dstPtr,scale,nProbe);
Jim_Unitary(prFilter,dstPtr,scale,nWork);
return TRUE;
}
/**************************************************************************************
I/Q相位校正算法 by Jim Fang at 2007
**************************************************************************************/
void WINAPI Jim_Pharev(void *pIvalue,void *pQvalue,void *pRev)
{
double *iPtr =(double *)pIvalue;
double *qPtr =(double *)pQvalue;
double *rPtr =(double *)pRev;
//计算自相关函数和互相关函数
double Rii=0.0,Rqq=0.0,Riq=0.0;
int i;
for ( i=0;i<fft_point;i++)
{
Rii=Rii+double(iPtr[i]*iPtr[i]);
Rqq=Rqq+double(qPtr[i]*qPtr[i]);
Riq=Riq+double(iPtr[i]*qPtr[i]);
}
rPtr[0]=sqrt(Rqq/Rii)-1;
rPtr[1]=asin(Riq/sqrt(Rii*Rqq));
}
/**************************************************************************************
能量归一化算法 by Jim Fang at 2007
FFT的源数据为WORD型，范围在±32768之间，用于显示能量的色阶范围为0-255
1.消除奇异点
2.找出峰值
3.归一化
**************************************************************************************/
void WINAPI Jim_Unitary(void *pSrcData,BYTE *pDstData,int nScale,int nWork)
{
//平滑处理后取最大值作为255，进行归一化。
int i,j;
double mod=0.0;
double max=0.0;
double sum=0.0;
int funny[10],nFunny;
double *srcPtr =(double *)pSrcData;
BYTE *dstPtr =(BYTE *)pDstData;
int nUnitary[9];
int nUnitary2mPw[9]={200,500,1000,1500,2000,2500,4000,6000,8000};
int nUnitary2mCw[9]={80,180,360,500,1000,1200,1500,1800,6000},
nUnitary4mCw[9];
//nUnitary2mPw[9]={200,500,1000,1500,2000,2500,4000,6000,8000};
//nUnitary4mCw[9]={80,180,360,500,1000,1200,1500,1800,6000};
//nProbe=2;
switch(nWork)
{
case 0:
for(i=0;i<9;i++)
{
nUnitary[i]=nUnitary2mPw[i];
}
break;
case 1:
for(i=0;i<9;i++)
{
nUnitary[i]=nUnitary2mCw[i];
}
break;
case 2:
for(i=0;i<9;i++)
{
nUnitary[i]=nUnitary2mCw[i];
}
break;
default:
for(i=0;i<9;i++)
{
nUnitary[i]=nUnitary2mPw[i];
}
break;
}
//奇异点处理，如果为奇异带如何处理
switch(nScale)
{
case 0:nFunny=5;funny[0]=0;funny[1]=10,funny[2]=51;funny[3]=205;funny[4]=246;break;
case 1:nFunny=5;funny[0]=0;funny[1]=41,funny[2]=82;funny[3]=174;funny[4]=215;break;
case 2:nFunny=10;funny[0]=0;funny[1]=72,funny[2]=86;funny[3]=91;funny[4]=99;
funny[5]=113;funny[6]=157;funny[7]=164;funny[8]=171;funny[9]=187;break;
case 3:nFunny=11;funny[0]=0;funny[1]=10,funny[2]=21;funny[3]=75;funny[4]=85;
funny[5]=96;funny[6]=150,funny[7]=171;funny[8]=181;funny[9]=236;
funny[10]=246;break;
default:nFunny=5;funny[0]=0;funny[1]=41,funny[2]=82;funny[3]=174;funny[4]=215;break;
}
for ( i=0;i<fft_point;i++ )
{
mod=srcPtr[i];
// 分段线性压缩，线性压缩是否最好？应采用二次函数压缩
if ( mod < nUnitary[0] )
mod -= nUnitary[0]/4;
if ( mod<0 ) mod = 0;
if ( mod < nUnitary[1] )
mod = mod/fft_divide;
else if ( mod < nUnitary[2] )
mod = mod*2.0/fft_divide;
else if ( mod < nUnitary[3] )
mod = mod*1.8/fft_divide;
else if ( mod < nUnitary[4] )
mod = mod*1.6/fft_divide;
else if ( mod < nUnitary[5] )
mod = mod*1.4/fft_divide;
else if ( mod < nUnitary[6] )
mod = mod*1.2/fft_divide;
else if ( mod < nUnitary[7] )
mod = mod*1.0/fft_divide;
else if ( mod < nUnitary[8] )
mod = mod*1.0/fft_divide;
else
mod = mod/(fft_divide);
if ( mod>=255 )
dstPtr[(i+fft_point/2)%fft_point] = 255;
else if ( mod<=0 )
dstPtr[(i+fft_point/2)%fft_point] = 0;
else
dstPtr[(i+fft_point/2)%fft_point] = (BYTE)((mod));
//滤除奇异点
for(j=0;j<nFunny;j++)
{
if(i==funny[j] )
dstPtr[(i+fft_point/2)%fft_point]=dstPtr[(i+fft_point/2)%fft_point-1];
}
}
//sum=sum/(256-nFunny);
//周期延拓
/*for ( i=0;i<fft_point;i++ )
{
if(i==funny[0] || i==funny[1]|| i==funny[2] || i==funny[3] || i==funny[4])
dstPtr[(i+fft_point/2)%fft_point]=dstPtr[(i+fft_point/2)%fft_point-1];
else
dstPtr[(i+fft_point/2)%fft_point] =srcPtr[i]*255.0/max;//sum*255.0/max;
}*/
}
/**************************************************************************************
点包络算法 by Jim Fang at 2007 (双向包络)
对256色的灰阶图象进行点包络，fft数据要进行奇异点滤除和平滑处理
所有点数始终为0-255来索引和输出，函数返回值为相对于基线的点数距离
返回值为16位数，高8位为正向包络，低8位为负向包络
先根据力矩法求质心位置，加速包络搜索速度
包络点即为：频谱的最大频率点(最好的算法为利用小波来计算)
offset:在屏幕坐标而言，上+下-
**************************************************************************************/
WORD WINAPI Jim_Cover(void *pSrcData,int base,int offset)
{
unsigned char peakPlus,peakNeg;
unsigned char startPlus,startNeg;
int cntPlus,cntNeg,i,j=0;
unsigned char *srcPtr = (unsigned char*)pSrcData;
bool bPlus,bNeg;
WORD cover;
//1.平滑，滤除噪声点，但不能消除真实的频谱信息
double smooth[256];
Jim_Smooth(256,srcPtr,smooth);
//2.阈值，利用力矩法求质心与 增益无关，取能量均值，考虑要取正向+反向的两个能量均值
double sumPlus=0.0,sumNeg=0.0;
double meanPlus,meanNeg;
double sumPowerPlus=0.0,sumPowerNeg=0.0;
int nCenterPlus,nCenterNeg;
for(i=0;i<256;i++)
{
if(i<base+offset)
{
sumNeg=sumNeg+smooth[i];
sumPowerNeg=sumPowerNeg+smooth[i]*(base+offset-i);
}
if(i>base+offset)
{
sumPlus=sumPlus+smooth[i];
sumPowerPlus=sumPowerPlus+smooth[i]*(i-(base+offset));
}
}
nCenterPlus=int(sumPowerPlus/sumPlus);
nCenterNeg=int(sumPowerNeg/sumNeg);
meanNeg=sumNeg/(base+offset);
meanPlus=sumPlus/(255-(base+offset));
double theshold=255*0.05;
//当阈值小于特定值后，不再包络
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -