📄 fftdoc.cpp
字号:
void CFFTDoc::Square()
{ InitArray();
for (i = 0; i < N/8; i++) {
samp[i].real=1;
samp[N/2+i].real=1;
samp[N/4+i].real=1;
samp[N*3/4+i].real=1;
}
CallFFT();
UpdateAllViews(NULL);
}
void CFFTDoc::CallFFT()
{ for (i = 0; i < N; i++)//幅度
tamp[i] = sqrt(samp[i].real*samp[i].real+samp[i].imag*samp[i].imag);
/* Display the spectrum curve计算最大值 */
tym = 1.0e-90;
for(i = 0; i < N; i++)
if (tamp[i] > tym) tym = tamp[i];
fft(samp,log(N)/log(2));
for (i = 0; i < N; i++)//幅度
famp[i] = sqrt(samp[i].real*samp[i].real+samp[i].imag*samp[i].imag);
/* Display the spectrum curve计算最大值 */
fym = 1.0e-90;
for(i = 0; i < N; i++)
if (famp[i] > fym) fym = famp[i];
}
//long CFFTDoc::OnMsgWave(WPARAM wParam, LPARAM lParam)
//{AfxMessageBox(_T("Doc"));
//return 0;
//}
void CFFTDoc::InitArray()
{ N = 512; /* FFT 点数 */
M = 1024; /* ZFFT 点数 */
for (i = 0; i < N; i++) {
samp[i].real = 0;
samp[i].imag = 0;
rsamp[i].real = 0;
rsamp[i].imag = 0;
tamp[i]=0;
famp[i]=0;
}
tym=0;fym=0;
}
void CFFTDoc::fSigFloat()
{ numav=16;estsize = (slice-ovlap)*(numav-1) + slice;m_fop=false;
InitArray();
double sig;
for (i=0; i<estsize; i++)
{sigfloat[i] = cos(6.0*1.0*PI*(double)i*2.0/slice);
sig=sigfloat[i];
}
for (i=0; i<estsize; i++)
{sigfloat[i] = sigfloat[i]+sin(35.0*1.0*PI*(double)i*2.0/slice);
sig=sigfloat[i];}
if(ntype==2) addnoise(sigfloat,estsize,0,0.2); /* 加入均匀噪声 */
else if(ntype==1) addnoise(sigfloat,estsize,1,0.2); /* 加入高斯噪声 */
if(m_bPower)
{CallPower(); }
else
{for (i = 0; i < N; i++) samp[i].real=sigfloat[i];
CallFFT();
}
UpdateAllViews(NULL);
}
/**************************************************************************
ADDNOISE - 在输入信号上加入均匀分布或高斯分布的噪声
输入参数:
double *Indata: 输入的未加噪声的信号序列;
unsigned length: 输入信号序列的长度;
int type : 加噪声的类型, 0:均匀 1:高斯
double nmult : 加入噪声的幅度值系数;
输出参数:
double *Indata: 在输入序列上加入噪声后,仍然将数据
放在 *Indata 数据区中.
*************************************************************************/
double * CFFTDoc::addnoise(double *Indata, unsigned int length, int type, double nmult)
{int i,j;
/* add noise to record */
if(type == 0)
for(i = 0 ; i < length ; i++) Indata[i] += nmult*funiform();
else
for(i = 0 ; i < length ; i++) Indata[i] += nmult*gaussian();
return(Indata);
}
/**************************************************************************
uniform - 产生零均值,均匀分布的随机数(-0.5 到 0.5 之间)
函数返回一个零均值,均匀分布的 double型随机数.
无参数, 调用 rand()函数.
double uniform()
*************************************************************************/
double CFFTDoc::funiform()
{//return 0;
return((double)(rand() & RAND_MAX) / RAND_MAX - 0.5);
}
/**************************************************************************
gaussian - 产生均值为零,单位方差的高斯分布随机数
函数返回 double型零均值,单位方差的高斯分布随机数.
使用 Box-Muller 转换方法,将一对均匀分布的随机变量映射为一对
高斯随机变量.
double gaussian()
*************************************************************************/
double CFFTDoc::gaussian()
{static int ready = 0; /* flag to indicated stored value */
static double gstore; /* place to store other value */
double v1,v2,r,fac,gaus;
/* make two numbers if none stored */
if(ready == 0) {
do {
v1 = 2.*funiform();
v2 = 2.*funiform();
r = v1*v1 + v2*v2;
} while(r > 1.0); /* make radius less than 1 */
/* remap v1 and v2 to two Gaussian numbers */
fac = sqrt(-2.*log(r)/r);
gstore = v1*fac; /* store one */
gaus = v2*fac; /* return one */
ready = 1; /* set ready flag */
}
else {
ready = 0; /* reset ready flag for next pair */
gaus = gstore; /* return the stored one */
}
return(gaus);
}
void CFFTDoc::CallPower()
{ int k,j,index,m;
double a,max;
double tempflt;
tym = 1.0e-90;
for(i = 0; i < estsize; i++)
{if ((sigfloat[i]>0)&&(sigfloat[i]>tym)) tym = sigfloat[i];
if ((sigfloat[i]<0)&&(-1.*sigfloat[i]>tym)) tym = -1.*sigfloat[i];
}
for (k=0; k<slice; k++) mag[k] = 0;
for (j=0; j<numav; j++){
for (k=0; k<slice; k++){
index = j*(slice-ovlap) + k;
samp[k].real = sigfloat[index];
samp[k].imag = 0;
}
switch (wtype) {
case 0: ham(samp,slice);
break;
case 1: han(samp,slice);
break;
case 2: triang(samp,slice);
break;
case 3: black(samp,slice);
break;
case 4: harris(samp,slice);
break;
default:
break;
}
m = log(slice)/log(2);
fft(samp,m);
a = (double) slice*slice;
for (k=0; k<slice; k++){
tempflt = samp[k].real * samp[k].real;
tempflt += samp[k].imag * samp[k].imag;
tempflt = tempflt/a;
mag[k] += tempflt;
}
}
/* Take log after averaging the magnitudes. */
for (k=0; k<slice/2; k++){
mag[k] = mag[k]/numav;
mag[k] = 10*log10(MAX(mag[k],1.e-14));
}
fym = 1.0e-90;
for(i = 0; i < slice/2; i++)
{if ((mag[i]>0)&&(mag[i]>fym)) fym = mag[i];
if ((mag[i]<0)&&(-1.*mag[i]>fym)) fym = -1.*mag[i];
}
j=0;
max=mag[0];
for(i = 1; i < slice/2; i++)
{if(mag[i]>=max){j=i;max=mag[i];}
}
m_cfreq=j*FRATE;
m_chigh=m_cfreq*HRATE;
}
/*************************************************************************
ham - Hamming 窗函数
输入参数:
COMPLEX *x : 输入复数数组指针;
int n : 输入数组长度.
输出参数:
加窗后结果放在输入数组中返回,
无输出参数.
void ham(COMPLEX *x, int n)
*************************************************************************/
void CFFTDoc::ham(COMPLEX *x, int n)
{
int i;
double ham,factor;
factor = 8.0*atan(1.0)/(n-1);
for (i = 0 ; i < n ; i++){
ham = 0.54 - 0.46*cos(factor*i);
x->real *= ham;
x->imag *= ham;
x++;
}
}
/*************************************************************************
han - Hanning 窗函数
输入参数:
COMPLEX *x : 输入复数数组指针;
int n : 输入数组长度.
输出参数:
加窗后结果放在输入数组中返回,
无输出参数.
void han(COMPLEX *x, int n)
*************************************************************************/
void CFFTDoc::han(COMPLEX *x, int n)
{int i;
double factor,han;
factor = 8.0*atan(1.0)/(n-1);
for (i = 0 ; i < n ; i++){
han = 0.5 - 0.5*cos(factor*i);
x->real *= han;
x->imag *= han;
x++;
}
}
/*************************************************************************
triang - triangle 窗函数
输入参数:
COMPLEX *x : 输入复数数组指针;
int n : 输入数组长度.
输出参数:
加窗后结果放在输入数组中返回,
无输出参数.
void triang(COMPLEX *,int n)
*************************************************************************/
void CFFTDoc::triang(COMPLEX *x, int n)
{
int i;
double tri,a;
a = 2.0/(n-1);
for (i = 0 ; i <= (n-1)/2 ; i++) {
tri = i*a;
x->real *= tri;
x->imag *= tri;
x++;
}
for ( ; i < n ; i++) {
tri = 2.0 - i*a;
x->real *= tri;
x->imag *= tri;
x++;
}
}
/*************************************************************************
black - Blackman 窗函数
输入参数:
COMPLEX *x : 输入复数数组指针;
int n : 输入数组长度.
输出参数:
加窗后结果放在输入数组中返回,
无输出参数.
void black(COMPLEX *x, int n)
*************************************************************************/
void CFFTDoc::black(COMPLEX *x, int n)
{int i;
double black,factor;
factor = 8.0*atan(1.0)/(n-1);
for (i=0; i<n; ++i){
black = 0.42 - 0.5*cos(factor*i) + 0.08*cos(2*factor*i);
x->real *= black;
x->imag *= black;
x++;
}
}
/*************************************************************************
harris - 4 term Blackman-Harris 窗函数
输入参数:
COMPLEX *x : 输入复数数组指针;
int n : 输入数组长度.
输出参数:
加窗后结果放在输入数组中返回,
无输出参数.
void harris(COMPLEX *x, int n)
*************************************************************************/
void CFFTDoc::harris(COMPLEX *x, int n)
{int i;
double harris,factor,arg;
factor = 8.0*atan(1.0)/n;
for (i=0; i<n; ++i){
arg = factor * i;
harris = 0.35875 - 0.48829*cos(arg) + 0.14128*cos(2*arg)
- 0.01168*cos(3*arg);
x->real *= harris;
x->imag *= harris;
x++;
}
}
void CFFTDoc::OnFileOpen()
{ numav=8;estsize = (slice-ovlap)*(numav-1) + slice;
InitArray();
if(m_fop)
{CallPower();
UpdateAllViews(NULL);
return;
}
m_fop=true;
// TODO: Add your command handler code here
CFileDialog *filedlg;
FILE *fp;
int index=0;
int h,m,s,hr;
float value;
filedlg = (CFileDialog*)new CFileDialog(TRUE,_T("out"),_T("Data"),NULL,"Date Save Files(*.out) | *.out|data Files (*.Dt)|*.Dt|All Files (*.*)|*.*||", NULL);
filedlg->DoModal();
CString lppathname=filedlg->GetPathName();
if((fp = fopen(lppathname, "rb" ))!=NULL)
{
//fscanf(fp,"\nhigh:%10d\n",&m_high);//高度
//while(m_high>=1)
// {
// fscanf(fp,"time:%2d:%2d:%2d.%2d\n",&h,&m,&s,&hr);
//UpdateData(false);
// ::Sleep(100);
for(index=0;index<NUM;index++)
{
//fscanf(fp,"%f \n",&value);
fscanf(fp,"%f,",&value);
sigfloat[index]=value;
//(*waveform)[index]=(double)value;
}
CallPower();
UpdateAllViews(NULL);
fscanf(fp,"\nhigh:%10d\n",&m_high);
//}
fclose(fp);
}
}
void CFFTDoc::Closefp()
{
m_fop=false;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -